####In Search of Self-censorship - R Analysis File####
###Xiaoxiao Shen, Ph.D. Candidate, Princeton University Department of Politics, xiaoxiao@princeton.edu###
###Rory Truex, Assistant Professor, Princeton University Department of Politics and Woodrow Wilson School of Public and International Affairs, rtruex@princeton.edu###

###LOAD PACKAGES AND FUNCTIONS###
setwd('/Users/rtruex/Google Drive/In Search of Preference Falsification/BJPOLS Revision')
rm(list=ls(all=TRUE))

library(countrycode)
library(ggplot2)
library(ggridges)
library(ggjoy)
library(survival)
library(Zelig)
library(foreign)
library(purrr)
library(stringr)
library(Rcpp)
library(Amelia)
library("GGally")
library(grid)
library("MissMech")
library("ggExtra")
library("gridExtra")


substrRight <- function(x, n){
  sapply(x, function(xx)
    substr(xx, (nchar(xx)-n+1), nchar(xx))
  )
}

require(grid)
vp.layout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)
arrange_ggplot2 <- function(..., nrow=NULL, ncol=NULL, as.table=FALSE) {
  dots <- list(...)
  n <- length(dots)
  if(is.null(nrow) & is.null(ncol)) { nrow = floor(n/2) ; ncol = ceiling(n/nrow)}
  if(is.null(nrow)) { nrow = ceiling(n/ncol)}
  if(is.null(ncol)) { ncol = ceiling(n/nrow)}
  
  grid.newpage()
  pushViewport(viewport(layout=grid.layout(nrow,ncol) ) )
  ii.p <- 1
  for(ii.row in seq(1, nrow)){
    ii.table.row <- ii.row  
    if(as.table) {ii.table.row <- nrow - ii.table.row + 1}
    for(ii.col in seq(1, ncol)){
      ii.table <- ii.p
      if(ii.p > n) break
      print(dots[[ii.table]], vp=vp.layout(ii.table.row, ii.col))
      ii.p <- ii.p + 1
    }
  }
}

###LOAD DATA###
wvs<-read.csv("In Search of Preference Falsification - WVS Longitudinal Data.csv",stringsAsFactors=FALSE)
data<-read.csv("data.wvs.full.csv",stringsAsFactors=FALSE)

###CREATE DATASET###

wvs$n<-1

wvs$conf.gov.obs<-1
wvs$conf.gov.obs[wvs$E069_11=="Not asked in survey"]<-0
wvs$conf.gov.dkna<-NA
wvs$conf.gov.dkna[wvs$E069_11=="Not asked in survey"]<-NA
wvs$conf.gov.dkna[wvs$E069_11=="Don´t know"]<-1
wvs$conf.gov.dkna[wvs$E069_11=="No answer"]<-1
wvs$conf.gov.dkna[wvs$E069_11=="A great deal"]<-0
wvs$conf.gov.dkna[wvs$E069_11=="None at all"]<-0
wvs$conf.gov.dkna[wvs$E069_11=="Not very much"]<-0
wvs$conf.gov.dkna[wvs$E069_11=="Quite a lot"]<-0
wvs$conf.gov<-NA
wvs$conf.gov[wvs$E069_11=="Not asked in survey"]<-NA
wvs$conf.gov[wvs$E069_11=="Don´t know"]<-NA
wvs$conf.gov[wvs$E069_11=="No answer"]<-NA 
wvs$conf.gov[wvs$E069_11=="A great deal"]<-4
wvs$conf.gov[wvs$E069_11=="None at all"]<-1
wvs$conf.gov[wvs$E069_11=="Not very much"]<-2
wvs$conf.gov[wvs$E069_11=="Quite a lot"]<-3
ftable(wvs$conf.gov.obs)
ftable(wvs$conf.gov.dkna)

wvs$resp.hr.obs<-1
wvs$resp.hr.obs[wvs$E124=="Not asked in survey"]<-0
wvs$resp.hr.dkna<-NA
wvs$resp.hr.dkna[wvs$E124=="Not asked in survey"]<-NA
wvs$resp.hr.dkna[wvs$E124=="Don´t know"]<-1 
wvs$resp.hr.dkna[wvs$E124=="No answer"]<-1
wvs$resp.hr.dkna[wvs$E124=="There is a lot of respect for individual human rights"]<-0
wvs$resp.hr.dkna[wvs$E124=="There is not much respect"]<-0
wvs$resp.hr.dkna[wvs$E124=="There is some respect"]<-0
wvs$resp.hr.dkna[wvs$E124=="There is no respect at all"]<-0
wvs$resp.hr<-NA
wvs$resp.hr[wvs$E124=="Not asked in survey"]<-NA
wvs$resp.hr[wvs$E124=="Don´t know"]<-NA
wvs$resp.hr[wvs$E124=="No answer"]<-NA 
wvs$resp.hr[wvs$E124=="There is a lot of respect for individual human rights"]<-4
wvs$resp.hr[wvs$E124=="There is not much respect"]<-2
wvs$resp.hr[wvs$E124=="There is some respect"]<-3
wvs$resp.hr[wvs$E124=="There is no respect at all"]<-1
ftable(wvs$resp.hr.obs)
ftable(wvs$resp.hr.dkna)

wvs$democraticness.obs<-1
wvs$democraticness.obs[wvs$E236=="Not asked in survey"]<-0
wvs$democraticness.dkna<-NA
wvs$democraticness.dkna[wvs$E236=="Not asked in survey"]<-NA
wvs$democraticness.dkna[wvs$E236=="Don´t know"]<-1
wvs$democraticness.dkna[wvs$E236=="No answer"]<-1
wvs$democraticness.dkna[wvs$E236=="Not at all democratic"]<-0
wvs$democraticness.dkna[wvs$E236=="2"]<-0
wvs$democraticness.dkna[wvs$E236=="3"]<-0
wvs$democraticness.dkna[wvs$E236=="4"]<-0
wvs$democraticness.dkna[wvs$E236=="5"]<-0
wvs$democraticness.dkna[wvs$E236=="6"]<-0
wvs$democraticness.dkna[wvs$E236=="7"]<-0
wvs$democraticness.dkna[wvs$E236=="8"]<-0
wvs$democraticness.dkna[wvs$E236=="9"]<-0
wvs$democraticness.dkna[wvs$E236=="Completely democratic"]<-0
wvs$democraticness<-NA
wvs$democraticness[wvs$E236=="Not asked in survey"]<-NA
wvs$democraticness[wvs$E236=="Don´t know"]<-NA 
wvs$democraticness[wvs$E236=="No answer"]<-NA 
wvs$democraticness[wvs$E236=="Not at all democratic"]<-1
wvs$democraticness[wvs$E236=="2"]<-2
wvs$democraticness[wvs$E236=="3"]<-3
wvs$democraticness[wvs$E236=="4"]<-4
wvs$democraticness[wvs$E236=="5"]<-5
wvs$democraticness[wvs$E236=="6"]<-6
wvs$democraticness[wvs$E236=="7"]<-7
wvs$democraticness[wvs$E236=="8"]<-8
wvs$democraticness[wvs$E236=="9"]<-9
wvs$democraticness[wvs$E236=="Completely democratic"]<-10
ftable(wvs$democraticness.obs)
ftable(wvs$democraticness.dkna)

wvs$lifesat.obs<-1
wvs$lifesat.obs[wvs$A170=="Not asked in survey"]<-0
wvs$lifesat.dkna<-NA
wvs$lifesat.dkna[wvs$A170=="Not asked in survey"]<-NA
wvs$lifesat.dkna[wvs$A170=="Don´t know"]<-1
wvs$lifesat.dkna[wvs$A170=="No answer"]<-1
wvs$lifesat.dkna[wvs$A170=="Dissatisfied"]<-0
wvs$lifesat.dkna[wvs$A170=="2"]<-0
wvs$lifesat.dkna[wvs$A170=="3"]<-0
wvs$lifesat.dkna[wvs$A170=="4"]<-0
wvs$lifesat.dkna[wvs$A170=="5"]<-0
wvs$lifesat.dkna[wvs$A170=="6"]<-0
wvs$lifesat.dkna[wvs$A170=="7"]<-0
wvs$lifesat.dkna[wvs$A170=="8"]<-0
wvs$lifesat.dkna[wvs$A170=="9"]<-0
wvs$lifesat.dkna[wvs$A170=="Satisfied"]<-0
ftable(wvs$lifesat.obs)
ftable(wvs$lifesat.dkna)

wvs$gentrust.obs<-1
wvs$gentrust.obs[wvs$A165=="Not asked in survey"]<-0
wvs$gentrust.dkna<-NA
wvs$gentrust.dkna[wvs$A165=="Not asked in survey"]<-NA
wvs$gentrust.dkna[wvs$A165=="Don´t know"]<-1
wvs$gentrust.dkna[wvs$A165=="No answer"]<-1
wvs$gentrust.dkna[wvs$A165=="Most people can be trusted"]<-0
wvs$gentrust.dkna[wvs$A165=="Can´t be too careful"]<-0
ftable(wvs$gentrust.obs)
ftable(wvs$gentrust.dkna)

wvs$conf.tv.obs<-1
wvs$conf.tv.obs[wvs$E069_10=="Not asked in survey"]<-0
wvs$conf.tv.dkna<-NA
wvs$conf.tv.dkna[wvs$E069_10=="Not asked in survey"]<-NA
wvs$conf.tv.dkna[wvs$E069_10=="Don´t know"]<-1
wvs$conf.tv.dkna[wvs$E069_10=="No answer"]<-1
wvs$conf.tv.dkna[wvs$E069_10=="A great deal"]<-0
wvs$conf.tv.dkna[wvs$E069_10=="None at all"]<-0
wvs$conf.tv.dkna[wvs$E069_10=="Not very much"]<-0
wvs$conf.tv.dkna[wvs$E069_10=="Quite a lot"]<-0
ftable(wvs$conf.tv.obs)
ftable(wvs$conf.tv.dkna)

vars<-c("E069_11","E124","E236","A170","A165","E069_10","A001", "A002", "A003", "A005", "A008", "A009", "B002", "B008", "C006", "D018", "D022", "D023", "D054", "D055", "D057", "D060", "E012")
  
wvs.agg<-aggregate(wvs$conf.gov.obs, by=list(wvs$S025), FUN='sum', na.rm=TRUE)[1]
for (j in vars) {
  try(eval(parse(text=paste("",sep=""))))
  try(eval(parse(text=paste("wvs$",j,".obs<-1",sep=""))))
  try(eval(parse(text=paste("wvs$",j,".obs[wvs$",j,"=='Not asked in survey']<-0",sep=""))))
  try(eval(parse(text=paste("wvs$",j,".dkna<-NA",sep=""))))
  try(eval(parse(text=paste("wvs$",j,".dkna[wvs$",j,"=='Not asked in survey']<-NA",sep=""))))
  try(eval(parse(text=paste("wvs$",j,".dkna[wvs$",j,"=='Don´t know']<-1",sep=""))))
  try(eval(parse(text=paste("wvs$",j,".dkna[wvs$",j,"=='No answer']<-1",sep=""))))
  try(eval(parse(text=paste("ftable(wvs$",j,".obs)",sep=""))))
  try(eval(parse(text=paste("ftable(wvs$",j,".dkna)",sep=""))))
  try(eval(parse(text=paste("wvs.agg",j,".ob<-aggregate(wvs$",j,".ob, by=list(wvs$S025), FUN='sum', na.rm=TRUE)[2]",sep=""))))
  try(eval(parse(text=paste("wvs.agg",j,".dkna<-aggregate(wvs$",j,".dkna, by=list(wvs$S025), FUN='sum', na.rm=TRUE)[2]",sep=""))))
  try(eval(parse(text=paste("colnames(wvs.agg",j,".ob)<-c('",j,".ob')",sep=""))))
  try(eval(parse(text=paste("colnames(wvs.agg",j,".dkna)<-c('",j,".dkna')",sep=""))))
  try(eval(parse(text=paste("wvs.agg<-cbind(wvs.agg, wvs.agg",j,".ob)",sep=""))))
  try(eval(parse(text=paste("wvs.agg<-cbind(wvs.agg, wvs.agg",j,".dkna)",sep=""))))
  try(eval(parse(text=paste("wvs.agg$",j,".dkna.rate<-(wvs.agg$",j,".dkna)/(wvs.agg$",j,".ob)",sep=""))))
  }

colnames(wvs.agg)[1]<-"country.year"
wvs.agg<-data.frame(wvs.agg)

data<-merge(data,wvs.agg,by="country.year",all.x=TRUE,all.y=FALSE)

###FULL DATA ANALYSIS###

##Imputation Model## - Get all figures working, then update this

data.impute<-data.frame(cbind(data$cow.year, data$country.year, data$country.name, data$year, data$n, data$polity, data$authoritarian, data$duration, data$executive.comp, data$party.comp, data$military, data$gdppc, data$oil,  data$repression, data$pts, data$edu, data$urbpop, data$E069_11.dkna.rate, data$E124.dkna.rate, data$E236.dkna.rate, data$A170.dkna.rate, data$A165.dkna.rate, data$E069_10.dkna.rate, data$A001.dkna.rate, data$A002.dkna.rate, data$A003.dkna.rate, data$A005.dkna.rate, data$A008.dkna.rate, data$A009.dkna.rate, data$B002.dkna.rate, data$B008.dkna.rate, data$C006.dkna.rate, data$D018.dkna.rate, data$D022.dkna.rate, data$D023.dkna.rate, data$D054.dkna.rate, data$D055.dkna.rate, data$D057.dkna.rate, data$D060.dkna.rate, data$E012.dkna.rate)) 
colnames(data.impute)<-c("cow.year", "country.year", "country.name", "year", "n", "polity","authoritarian", "duration", "executive.comp", "party.comp", "military", "gdppc", "oil", "repression", "pts", "edu", "urbpop", "E069_11.dkna.rate", "E124.dkna.rate", "E236.dkna.rate", "A170.dkna.rate", "A165.dkna.rate", "E069_10.dkna.rate", "A001.dkna.rate", "A002.dkna.rate", "A003.dkna.rate", "A005.dkna.rate", "A008.dkna.rate", "A009.dkna.rate", "B002.dkna.rate", "B008.dkna.rate", "C006.dkna.rate", "D018.dkna.rate", "D022.dkna.rate", "D023.dkna.rate", "D054.dkna.rate", "D055.dkna.rate", "D057.dkna.rate", "D060.dkna.rate", "E012.dkna.rate")

data.impute$n<-as.numeric(paste(data.impute$n))
data.impute$polity<-as.numeric(paste(data.impute$polity))           
data.impute$gdppc<-as.numeric(paste(data.impute$gdppc))           
data.impute$pts<-as.numeric(paste(data.impute$pts))
data.impute$edu<-as.numeric(paste(data.impute$edu))
data.impute$E069_11.dkna.rate<-as.numeric(paste(data.impute$E069_11.dkna.rate))
data.impute$E124.dkna.rate<-as.numeric(paste(data.impute$E124.dkna.rate))
data.impute$E236.dkna.rate<-as.numeric(paste(data.impute$E236.dkna.rate))
data.impute$A170.dkna.rate<-as.numeric(paste(data.impute$A170.dkna.rate))
data.impute$A165.dkna.rate<-as.numeric(paste(data.impute$A165.dkna.rate))
data.impute$E069_10.dkna.rate<-as.numeric(paste(data.impute$E069_10.dkna.rate))
data.impute$A001.dkna.rate<-as.numeric(paste(data.impute$A001.dkna.rate))
data.impute$A002.dkna.rate<-as.numeric(paste(data.impute$A002.dkna.rate))   
data.impute$A003.dkna.rate<-as.numeric(paste(data.impute$A003.dkna.rate))
data.impute$A005.dkna.rate<-as.numeric(paste(data.impute$A005.dkna.rate))
data.impute$A008.dkna.rate<-as.numeric(paste(data.impute$A008.dkna.rate))
data.impute$A009.dkna.rate<-as.numeric(paste(data.impute$A009.dkna.rate))
data.impute$B002.dkna.rate<-as.numeric(paste(data.impute$B002.dkna.rate))
data.impute$B008.dkna.rate<-as.numeric(paste(data.impute$B008.dkna.rate))   
data.impute$C006.dkna.rate<-as.numeric(paste(data.impute$C006.dkna.rate))
data.impute$D018.dkna.rate<-as.numeric(paste(data.impute$D018.dkna.rate))
data.impute$D022.dkna.rate<-as.numeric(paste(data.impute$D022.dkna.rate))
data.impute$D023.dkna.rate<-as.numeric(paste(data.impute$D023.dkna.rate))
data.impute$D054.dkna.rate<-as.numeric(paste(data.impute$D054.dkna.rate))
data.impute$D055.dkna.rate<-as.numeric(paste(data.impute$D055.dkna.rate))   
data.impute$D057.dkna.rate<-as.numeric(paste(data.impute$D057.dkna.rate))
data.impute$D060.dkna.rate<-as.numeric(paste(data.impute$D060.dkna.rate))
data.impute$E012.dkna.rate<-as.numeric(paste(data.impute$E012.dkna.rate))
data.impute$duration<-as.numeric(paste(data.impute$duration))
data.impute$gdppc<-as.numeric(paste(data.impute$gdppc))
data.impute$repression<-as.numeric(paste(data.impute$repression))
data.impute$edu<-as.numeric(paste(data.impute$edu))
data.impute$urbpop<-as.numeric(paste(data.impute$urbpop))
data.impute$party.comp<-as.numeric(paste(data.impute$party.comp))
data.impute$military<-as.numeric(paste(data.impute$military))
data.impute$executive.comp<-as.numeric(paste(data.impute$executive.comp))
data.impute$oil<-as.numeric(paste(data.impute$oil))
data.impute$pts<-as.numeric(paste(data.impute$pts))

data.impute$exclude<-0
data.impute$exclude[(is.na(data.impute$E069_11.dkna.rate)=="TRUE")&(is.na(data.impute$E124.dkna.rate)=="TRUE")&(is.na(data.impute$E236.dkna.rate)=="TRUE")]<-1
data.impute$exclude[(data.impute$E069_11.dkna.rate==0)&(data.impute$E124.dkna.rate==0)&(data.impute$E236.dkna.rate==0)&(data.impute$A170.dkna.rate==0)&(data.impute$A165.dkna.rate==0)&(data.impute$E069_10.dkna.rate==0)]<-1
data.impute$exclude[(is.na(data.impute$polity)=="TRUE")]<-1
data.impute<-subset(data.impute,data.impute$exclude==0)

#Testing MCAR

data.impute.mcartest<-data.impute[,c(18,19,20,21,22,23)] #limiting test to 6 core missing variables
TestMCARNormality(data = data.impute.mcartest, del.lesscases = 1)

#Imputation

a.out <- amelia(data.impute, m = 50,idvars = c("cow.year","country.year","country.name","year","n","polity","duration","executive.comp","party.comp","military","gdppc","oil","repression","pts","edu","urbpop","exclude"), noms="authoritarian", emburn=c(1,100000))

a.out <- transform(a.out,preffalse.ind=((E069_11.dkna.rate + E124.dkna.rate + E236.dkna.rate)/3)-((E069_10.dkna.rate + A165.dkna.rate + A170.dkna.rate)/3))
a.out <- transform(a.out,preffalse.ind.min=E069_11.dkna.rate-A165.dkna.rate)
a.out <- transform(a.out,preffalse.ind.mid=(((E069_11.dkna.rate + E236.dkna.rate)/2)-((A165.dkna.rate + A170.dkna.rate)/2)))
a.out <- transform(a.out,preffalse.unadj.ind=((E069_11.dkna.rate + E124.dkna.rate + E236.dkna.rate)/3))
a.out <- transform(a.out,preffalse.unadj.ind.min=E069_11.dkna.rate)
a.out <- transform(a.out,preffalse.unadj.ind.mid=(((E069_11.dkna.rate + E236.dkna.rate)/2)))
a.out <- transform(a.out,preffalse.ind.20=((E069_11.dkna.rate + E124.dkna.rate + E236.dkna.rate)/3)-((E069_10.dkna.rate + A165.dkna.rate + A170.dkna.rate + A001.dkna.rate + A002.dkna.rate + A003.dkna.rate + A005.dkna.rate + A008.dkna.rate + A009.dkna.rate + B002.dkna.rate + B008.dkna.rate + C006.dkna.rate + D018.dkna.rate + D022.dkna.rate + D023.dkna.rate + D054.dkna.rate + D055.dkna.rate + D057.dkna.rate + D060.dkna.rate + E012.dkna.rate)/20))

summary(a.out)

write.amelia(obj=a.out, file.stem = "outdata")

for (j in 1:50) {
  try(eval(parse(text=paste("",sep=""))))
  try(eval(parse(text=paste("outdata",j,"<-read.csv('outdata",j,".csv')",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$sens.ob<-3*outdata",j,"$n",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$sens.dkna.rate<-(outdata",j,"$E069_11.dkna.rate + outdata",j,"$E124.dkna.rate + outdata",j,"$E236.dkna.rate)/3",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$sens.dkna.se<-sqrt((abs(outdata",j,"$sens.dkna.rate*(1-outdata",j,"$sens.dkna.rate)))/outdata",j,"$sens.ob)",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$sens.dkna.u95ci<-outdata",j,"$sens.dkna.rate+1.96*outdata",j,"$sens.dkna.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$sens.dkna.l95ci<-outdata",j,"$sens.dkna.rate-1.96*outdata",j,"$sens.dkna.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.ob<-3*outdata",j,"$n",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.dkna.rate<-(outdata",j,"$E069_10.dkna.rate + outdata",j,"$A165.dkna.rate + outdata",j,"$A170.dkna.rate)/3",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.dkna.se<-sqrt((abs(outdata",j,"$nonsens.dkna.rate*(1-outdata",j,"$nonsens.dkna.rate)))/outdata",j,"$nonsens.ob)",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.dkna.u95ci<-outdata",j,"$nonsens.dkna.rate+1.96*outdata",j,"$nonsens.dkna.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.dkna.l95ci<-outdata",j,"$nonsens.dkna.rate-1.96*outdata",j,"$nonsens.dkna.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.20.ob<-20*outdata",j,"$n",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.20.dkna.rate<-(outdata",j,"$E069_10.dkna.rate + outdata",j,"$A165.dkna.rate + outdata",j,"$A170.dkna.rate + outdata",j,"$A001.dkna.rate + outdata",j,"$A002.dkna.rate + outdata",j,"$A003.dkna.rate + outdata",j,"$A005.dkna.rate + outdata",j,"$A008.dkna.rate + outdata",j,"$A009.dkna.rate + outdata",j,"$B002.dkna.rate + outdata",j,"$B008.dkna.rate + outdata",j,"$C006.dkna.rate + outdata",j,"$D018.dkna.rate + outdata",j,"$D022.dkna.rate + outdata",j,"$D023.dkna.rate + outdata",j,"$D054.dkna.rate + outdata",j,"$D055.dkna.rate + outdata",j,"$D057.dkna.rate + outdata",j,"$D060.dkna.rate + outdata",j,"$E012.dkna.rate)/20",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.20.dkna.se<-sqrt((abs(outdata",j,"$nonsens.20.dkna.rate*(1-outdata",j,"$nonsens.20.dkna.rate)))/outdata",j,"$nonsens.20.ob)",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.20.dkna.u95ci<-outdata",j,"$nonsens.20.dkna.rate+1.96*outdata",j,"$nonsens.20.dkna.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.20.dkna.l95ci<-outdata",j,"$nonsens.20.dkna.rate-1.96*outdata",j,"$nonsens.20.dkna.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind<-outdata",j,"$sens.dkna.rate-outdata",j,"$nonsens.dkna.rate",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.se<-sqrt((outdata",j,"$nonsens.dkna.se^2)+(outdata",j,"$sens.dkna.se^2))",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.u95ci<-outdata",j,"$preffalse.ind+1.96*outdata",j,"$preffalse.ind.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.l95ci<-outdata",j,"$preffalse.ind-1.96*outdata",j,"$preffalse.ind.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.min<-outdata",j,"$E069_11.dkna.rate-outdata",j,"$E069_10.dkna.rate",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$sens.min.dkna.se<-sqrt((abs(outdata",j,"$E069_11.dkna.rate*(1-outdata",j,"$E069_11.dkna.rate)))/outdata",j,"$n)",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$nonsens.min.dkna.se<-sqrt((abs(outdata",j,"$E069_10.dkna.rate*(1-outdata",j,"$E069_10.dkna.rate)))/outdata",j,"$n)",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.min.se<-sqrt((outdata",j,"$nonsens.min.dkna.se^2)+(outdata",j,"$sens.min.dkna.se^2))",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.min.u95ci<-outdata",j,"$preffalse.ind.min+1.96*outdata",j,"$preffalse.ind.min.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.min.l95ci<-outdata",j,"$preffalse.ind.min-1.96*outdata",j,"$preffalse.ind.min.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.20<-outdata",j,"$sens.dkna.rate-outdata",j,"$nonsens.20.dkna.rate",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.20.se<-sqrt((outdata",j,"$nonsens.20.dkna.se^2)+(outdata",j,"$sens.dkna.se^2))",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.20.u95ci<-outdata",j,"$preffalse.ind.20+1.96*outdata",j,"$preffalse.ind.20.se",sep=""))))
  try(eval(parse(text=paste("outdata",j,"$preffalse.ind.20.l95ci<-outdata",j,"$preffalse.ind.20-1.96*outdata",j,"$preffalse.ind.20.se",sep=""))))
  }  

##Calculate Falsification Index##

outdata.agg<-rbind(outdata1, outdata2, outdata3, outdata4, outdata5, outdata6, outdata7, outdata8, outdata9, outdata10, outdata11, outdata12, outdata13, outdata14, outdata15, outdata16, outdata17, outdata18, outdata19, outdata20, outdata21, outdata22, outdata23, outdata24, outdata25, outdata26, outdata27, outdata28, outdata29, outdata30, outdata31, outdata32, outdata33, outdata34, outdata35, outdata36, outdata37, outdata38, outdata39, outdata40, outdata41, outdata42, outdata43, outdata44, outdata45, outdata46, outdata47, outdata48, outdata49, outdata50)

outdata.preffalse.ind<-aggregate(outdata.agg$preffalse.ind, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.ind)<-c("country.year","cow.year","country.name","preffalse.ind")

outdata.preffalse.ind.mid<-aggregate(outdata.agg$preffalse.ind.mid, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.ind.mid)<-c("country.year","cow.year","country.name","preffalse.ind.mid")

outdata.preffalse.ind.min<-aggregate(outdata.agg$preffalse.ind.min, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.ind.min)<-c("country.year","cow.year","country.name","preffalse.ind.min")

outdata.preffalse.ind.20<-aggregate(outdata.agg$preffalse.ind.20, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.ind.20)<-c("country.year","cow.year","country.name","preffalse.ind.20")

outdata.preffalse.unadj.ind<-aggregate(outdata.agg$preffalse.unadj.ind, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.unadj.ind)<-c("country.year","cow.year","country.name","preffalse.unadj.ind")

outdata.preffalse.unadj.ind.mid<-aggregate(outdata.agg$preffalse.unadj.ind.mid, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.unadj.ind.mid)<-c("country.year","cow.year","country.name","preffalse.unadj.ind.mid")

outdata.preffalse.unadj.ind.min<-aggregate(outdata.agg$preffalse.unadj.ind.min, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.unadj.ind.min)<-c("country.year","cow.year","country.name","preffalse.unadj.ind.min")

outdata.preffalse.ind.u95ci<-aggregate(outdata.agg$preffalse.ind.u95ci, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE) 
colnames(outdata.preffalse.ind.u95ci)<-c("country.year","cow.year","country.name","preffalse.ind.u95ci") 

outdata.preffalse.ind.l95ci<-aggregate(outdata.agg$preffalse.ind.l95ci, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.ind.l95ci)<-c("country.year","cow.year","country.name","preffalse.ind.l95ci")

outdata.preffalse.ind.20.u95ci<-aggregate(outdata.agg$preffalse.ind.20.u95ci, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE) 
colnames(outdata.preffalse.ind.20.u95ci)<-c("country.year","cow.year","country.name","preffalse.ind.20.u95ci") 

outdata.preffalse.ind.20.l95ci<-aggregate(outdata.agg$preffalse.ind.20.l95ci, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.preffalse.ind.20.l95ci)<-c("country.year","cow.year","country.name","preffalse.ind.20.l95ci")

outdata.nonsens.dkna.rate<-aggregate(outdata.agg$nonsens.dkna.rate, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.nonsens.dkna.rate)<-c("country.year","cow.year","country.name","nonsens.dkna.rate")

outdata.sens.dkna.rate<-aggregate(outdata.agg$sens.dkna.rate, by=list(outdata.agg$country.year,outdata.agg$cow.year,outdata.agg$country.name), FUN="mean", na.rm=TRUE)
colnames(outdata.sens.dkna.rate)<-c("country.year","cow.year","country.name","sens.dkna.rate")

outdata<-merge(outdata.preffalse.ind,outdata.preffalse.ind.u95ci,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.ind.l95ci,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.ind.20,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.ind.20.l95ci,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.ind.20.u95ci,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.ind.mid,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.ind.min,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.unadj.ind,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.unadj.ind.mid,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.preffalse.unadj.ind.min,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.sens.dkna.rate,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)
outdata<-merge(outdata,outdata.nonsens.dkna.rate,by=c("country.year", "cow.year", "country.name"), all.x=TRUE,all.y=FALSE)

##Figure - Preference Falsification Index and Regime Type##

data.impute<-merge(data.impute,outdata,by="cow.year",all.x=TRUE,all.y=FALSE) #Fix .x and .y problem
data.impute$regime<-NA
data.impute$regime[(data.impute$polity<5)==TRUE]<-"Authoritarian"
data.impute$regime[(data.impute$polity==10)==TRUE]<-"Full Democracy"
data.impute.comp<-subset(data.impute, is.na(data.impute$regime)==FALSE)

pdf('fig-density.pdf', width=6, height=4)                                                                                                                                                                                                                                                                                                                                                                                                                                              
ggplot() + geom_line(aes(x=preffalse.ind,colour = regime), trim=TRUE, adjust=2.5, size=1, stat="density",alpha=.6, data.impute.comp)  + guides(color=guide_legend(title=NULL)) + theme_classic() + xlim(c(-.05,.17))+ theme(axis.title.x = element_text(size = 12, vjust = .25)) + xlab("Self-censorship Index") + ylab("Density") + geom_vline(xintercept = 0, linetype="dashed", size=.3) + scale_color_manual(values=c("dodgerblue3","red1")) + theme(legend.position=c(0.8, 0.875)) + geom_rug(aes(x=preffalse.ind,colour = regime), alpha=.6, data.impute.comp) +  annotate("text", x = 0.1168591997, y = 4, label = "China (2007)", size=2.6) + geom_segment(aes(x = 0.1168591997, y = 0, xend = 0.1168591997, yend = 3.5), linetype="dotted", color="grey50") +  annotate("text", x = 0.1000000000, y = 6, label = "Morocco (2007)", size=2.6) +  geom_segment(aes(x = 0.1000000000, y = 0, xend = 0.1000000000, yend = 5.5), linetype="dotted", color="grey50")  +  annotate("text", x = 0.0940835812, y = 8, label = "Pakistan (2001)", size=2.6) + geom_segment(aes(x = 0.0940835812, y = 0, xend = 0.0940835812, yend = 7.5), linetype="dotted", color="grey50") + annotate("text", x = 0.0830612692, y = 10, label = "China (2001)", size=2.6) + geom_segment(aes(x = 0.0830612692, y = 0, xend = 0.0830612692, yend = 9.5), linetype="dotted", color="grey50") + annotate("text", x = 0.0750724638, y = 12, label = "China (2012)", size=2.6) + geom_segment(aes(x = 0.0750724638, y = 0, xend = 0.0750724638, yend = 11.5), linetype="dotted", color="grey50")
dev.off()

pdf('fig-density-20.pdf', width=6, height=4)                                                                                                                                                                                                                                                                                                                                                                                                                                              
ggplot() + geom_line(aes(x=preffalse.ind.20,colour = regime), trim=TRUE, adjust=2.5, size=1, stat="density",alpha=.6, data.impute.comp)  + guides(color=guide_legend(title=NULL)) + theme_classic() + xlim(c(-.05,.17))+ theme(axis.title.x = element_text(size = 12, vjust = .25)) + xlab("Self-censorship Index") + ylab("Density") + geom_vline(xintercept = 0, linetype="dashed", size=.3) + scale_color_manual(values=c("dodgerblue3","red1")) + theme(legend.position=c(0.8, 0.875)) + geom_rug(aes(x=preffalse.ind.20,colour = regime), alpha=.6, data.impute.comp) +  annotate("text", x = 0.0924075004, y = 8, label = "China (2007)", size=2.6) + geom_segment(aes(x = 0.0924075004, y = 0, xend = 0.0924075004, yend = 7.5), linetype="dotted", color="grey50") +  annotate("text", x = 0.0985277778, y = 6, label = "Morocco (2007)", size=2.6) +  geom_segment(aes(x = 0.0985277778, y = 0, xend = 0.0985277778, yend = 5.5), linetype="dotted", color="grey50")  +  annotate("text", x = 0.1698344361, y = 4, label = "Pakistan (2001)", size=2.6) + geom_segment(aes(x = 0.1698344361, y = 0, xend = 0.1698344361, yend = 3.5), linetype="dotted", color="grey50") + annotate("text", x = 0.0616112692, y = 14, label = "China (2001)", size=2.6) + geom_segment(aes(x = 0.0616112692, y = 0, xend = 0.0616112692, yend = 13.5), linetype="dotted", color="grey50") + annotate("text", x = 0.0657260164, y = 12, label = "China (2012)", size=2.6) + geom_segment(aes(x = 0.0657260164, y = 0, xend = 0.0657260164, yend = 11.5), linetype="dotted", color="grey50") + annotate("text", x = 0.0690134724, y = 10, label = "Burkina Faso (2007)", size=2.6) + geom_segment(aes(x = 0.0690134724, y = 0, xend = 0.0690134724, yend = 9.5), linetype="dotted", color="grey50")
dev.off()

mean(data.impute.comp$preffalse.ind)
sd(data.impute.comp$preffalse.ind)

mean(data.impute.comp$preffalse.ind[data.impute.comp$regime=="Full Democracy"])
sd(data.impute.comp$preffalse.ind[data.impute.comp$regime=="Full Democracy"])

mean(data.impute.comp$preffalse.ind[data.impute.comp$regime=="Authoritarian"])
sd(data.impute.comp$preffalse.ind[data.impute.comp$regime=="Authoritarian"])

##Table - Testing for Differences in Falsification##

m1<-zelig(preffalse.ind~polity, data=a.out, model="ls") 
summary(m1)
m2<-zelig(preffalse.unadj.ind~polity, data=a.out, model="ls") 
summary(m2)
m3<-zelig(preffalse.ind.20~polity, data=a.out, model="ls") 
summary(m3)

m4<-zelig(preffalse.ind~authoritarian, data=a.out, model="ls") 
summary(m4)
m5<-zelig(preffalse.unadj.ind~authoritarian, data=a.out, model="ls") 
summary(m5)
m6<-zelig(preffalse.ind.20~authoritarian, data=a.out, model="ls") 
summary(m6)

###AUTHORITARIAN REGIMES ANALYSIS###

data.impute<-subset(data.impute, data.impute$authoritarian==1)
data.impute$label<-paste(data.impute$country.name.x,data.impute$year, sep=" - " )

pdf('fig2-wvs-pfi-dotplot-upd.pdf', width=7.5, height=10.5)
p1<-ggplot(data.impute, aes(y=reorder(label,preffalse.ind), x=preffalse.ind)) + xlim(c(-.13,.13)) + xlab("Self-censorship Index") + ylab("Country-Year") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue3", size=2.5) + theme_classic() + geom_segment(aes(y = reorder(label,preffalse.ind), x = preffalse.ind.l95ci, xend = preffalse.ind.u95ci, yend = reorder(label,preffalse.ind)), color="dodgerblue3", alpha=.4, lwd=.8) + geom_vline(xintercept = 0.0,  lty="dashed", alpha=.6) #+ geom_rug(outside=TRUE)# + coord_cartesian(clip = "off")
p2<-ggplot(data.impute.comp, aes(y=preffalse.ind,x = regime, color=regime)) + geom_boxplot(alpha=.5) + coord_flip() + theme_classic() + theme(legend.position="none") + ylim(c(-.155,.13)) + xlab("Regime") + scale_color_manual(values=c("dodgerblue3","red1"), labels=c("Authoritarian","Full Democracy"))  +theme(plot.background = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), panel.border = element_blank(),panel.background = element_blank(),axis.title.x = element_blank(),axis.text.x = element_blank(),axis.ticks = element_blank(),axis.line = element_blank()) #
grid.arrange(p1, p2, ncol=1, nrow=2, heights=c(10, 1))
dev.off() 

pdf('fig2-wvs-pfi-dotplot-upd-20.pdf', width=7.5, height=10.5)
p1<-ggplot(data.impute, aes(y=reorder(label,preffalse.ind.20), x=preffalse.ind.20)) + xlim(c(-.13,.13)) + xlab("Self-censorship Index") + ylab("Country-Year") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue3", size=2.5) + theme_classic() + geom_segment(aes(y = reorder(label,preffalse.ind.20), x = preffalse.ind.20.l95ci, xend = preffalse.ind.20.u95ci, yend = reorder(label,preffalse.ind.20)), color="dodgerblue3", alpha=.4, lwd=.8) + geom_vline(xintercept = 0.0,  lty="dashed", alpha=.6) #+ geom_rug(outside=TRUE)# + coord_cartesian(clip = "off")
p2<-ggplot(data.impute.comp, aes(y=preffalse.ind.20,x = regime, color=regime)) + geom_boxplot(alpha=.5) + coord_flip() + theme_classic() + theme(legend.position="none") + ylim(c(-.155,.13)) + xlab("Regime") + scale_color_manual(values=c("dodgerblue3","red1"), labels=c("Authoritarian","Full Democracy"))  +theme(plot.background = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), panel.border = element_blank(),panel.background = element_blank(),axis.title.x = element_blank(),axis.text.x = element_blank(),axis.ticks = element_blank(),axis.line = element_blank()) #
grid.arrange(p1, p2, ncol=1, nrow=2, heights=c(10, 1))
dev.off() 

##Sensitivity Analysis##

m1<-zelig(preffalse.ind~executive.comp, data=a.out, model="ls") 
m2<-zelig(preffalse.ind~executive.comp + military + party.comp, data=a.out, model="ls") 
m3<-zelig(preffalse.ind~executive.comp + military + party.comp + repression + duration, data=a.out, model="ls") 
m4<-zelig(preffalse.ind~executive.comp + military + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m5<-zelig(preffalse.ind~military, data=a.out, model="ls") 
m6<-zelig(preffalse.ind~military+ executive.comp + party.comp, data=a.out, model="ls") 
m7<-zelig(preffalse.ind~military+ executive.comp  + party.comp + repression + duration, data=a.out, model="ls") 
m8<-zelig(preffalse.ind~military+ executive.comp + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m9<-zelig(preffalse.ind~party.comp, data=a.out, model="ls") 
m10<-zelig(preffalse.ind~party.comp+ executive.comp + military, data=a.out, model="ls") 
m11<-zelig(preffalse.ind~party.comp+ executive.comp  +military + repression + duration, data=a.out, model="ls") 
m12<-zelig(preffalse.ind~party.comp+ executive.comp + military + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m13<-zelig(preffalse.ind~repression, data=a.out, model="ls") 
m14<-zelig(preffalse.ind~repression+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m15<-zelig(preffalse.ind~repression + party.comp+ executive.comp  +military +duration, data=a.out, model="ls") 
m16<-zelig(preffalse.ind~repression+ party.comp+ executive.comp + military  + duration + oil + gdppc + edu, data=a.out, model="ls") 
m17<-zelig(preffalse.ind~duration, data=a.out, model="ls") 
m18<-zelig(preffalse.ind~duration+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m19<-zelig(preffalse.ind~duration+ party.comp+ executive.comp  +military + repression, data=a.out, model="ls") 
m20<-zelig(preffalse.ind~duration +party.comp+ executive.comp + military + repression + oil + gdppc + edu, data=a.out, model="ls") 
m21<-zelig(preffalse.ind.min~executive.comp, data=a.out, model="ls") 
m22<-zelig(preffalse.ind.min~executive.comp + military + party.comp, data=a.out, model="ls") 
m23<-zelig(preffalse.ind.min~executive.comp + military + party.comp + repression + duration, data=a.out, model="ls") 
m24<-zelig(preffalse.ind.min~executive.comp + military + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m25<-zelig(preffalse.ind.min~military, data=a.out, model="ls") 
m26<-zelig(preffalse.ind.min~military+ executive.comp + party.comp, data=a.out, model="ls") 
m27<-zelig(preffalse.ind.min~military+ executive.comp  + party.comp + repression + duration, data=a.out, model="ls") 
m28<-zelig(preffalse.ind.min~military+ executive.comp + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m29<-zelig(preffalse.ind.min~party.comp, data=a.out, model="ls") 
m30<-zelig(preffalse.ind.min~party.comp+ executive.comp + military, data=a.out, model="ls") 
m31<-zelig(preffalse.ind.min~party.comp+ executive.comp  +military + repression + duration, data=a.out, model="ls") 
m32<-zelig(preffalse.ind.min~party.comp+ executive.comp + military + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m33<-zelig(preffalse.ind.min~repression, data=a.out, model="ls") 
m34<-zelig(preffalse.ind.min~repression+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m35<-zelig(preffalse.ind.min~repression + party.comp+ executive.comp  +military +duration, data=a.out, model="ls") 
m36<-zelig(preffalse.ind.min~repression+ party.comp+ executive.comp + military  + duration + oil + gdppc + edu, data=a.out, model="ls") 
m37<-zelig(preffalse.ind.min~duration, data=a.out, model="ls") 
m38<-zelig(preffalse.ind.min~duration+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m39<-zelig(preffalse.ind.min~duration+ party.comp+ executive.comp  +military + repression, data=a.out, model="ls") 
m40<-zelig(preffalse.ind.min~duration +party.comp+ executive.comp + military + repression + oil + gdppc + edu, data=a.out, model="ls") 
m41<-zelig(preffalse.ind.mid~executive.comp, data=a.out, model="ls") 
m42<-zelig(preffalse.ind.mid~executive.comp + military + party.comp, data=a.out, model="ls") 
m43<-zelig(preffalse.ind.mid~executive.comp + military + party.comp + repression + duration, data=a.out, model="ls") 
m44<-zelig(preffalse.ind.mid~executive.comp + military + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m45<-zelig(preffalse.ind.mid~military, data=a.out, model="ls") 
m46<-zelig(preffalse.ind.mid~military+ executive.comp + party.comp, data=a.out, model="ls") 
m47<-zelig(preffalse.ind.mid~military+ executive.comp  + party.comp + repression + duration, data=a.out, model="ls") 
m48<-zelig(preffalse.ind.mid~military+ executive.comp + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m49<-zelig(preffalse.ind.mid~party.comp, data=a.out, model="ls") 
m50<-zelig(preffalse.ind.mid~party.comp+ executive.comp + military, data=a.out, model="ls") 
m51<-zelig(preffalse.ind.mid~party.comp+ executive.comp  +military + repression + duration, data=a.out, model="ls") 
m52<-zelig(preffalse.ind.mid~party.comp+ executive.comp + military + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m53<-zelig(preffalse.ind.mid~repression, data=a.out, model="ls") 
m54<-zelig(preffalse.ind.mid~repression+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m55<-zelig(preffalse.ind.mid~repression + party.comp+ executive.comp  +military +duration, data=a.out, model="ls") 
m56<-zelig(preffalse.ind.mid~repression+ party.comp+ executive.comp + military  + duration + oil + gdppc + edu, data=a.out, model="ls") 
m57<-zelig(preffalse.ind.mid~duration, data=a.out, model="ls") 
m58<-zelig(preffalse.ind.mid~duration+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m59<-zelig(preffalse.ind.mid~duration+ party.comp+ executive.comp  +military + repression, data=a.out, model="ls") 
m60<-zelig(preffalse.ind.mid~duration +party.comp+ executive.comp + military + repression + oil + gdppc + edu, data=a.out, model="ls") 
m61<-zelig(preffalse.unadj.ind~executive.comp, data=a.out, model="ls") 
m62<-zelig(preffalse.unadj.ind~executive.comp + military + party.comp, data=a.out, model="ls") 
m63<-zelig(preffalse.unadj.ind~executive.comp + military + party.comp + repression + duration, data=a.out, model="ls") 
m64<-zelig(preffalse.unadj.ind~executive.comp + military + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m65<-zelig(preffalse.unadj.ind~military, data=a.out, model="ls") 
m66<-zelig(preffalse.unadj.ind~military+ executive.comp + party.comp, data=a.out, model="ls") 
m67<-zelig(preffalse.unadj.ind~military+ executive.comp  + party.comp + repression + duration, data=a.out, model="ls") 
m68<-zelig(preffalse.unadj.ind~military+ executive.comp + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m69<-zelig(preffalse.unadj.ind~party.comp, data=a.out, model="ls") 
m70<-zelig(preffalse.unadj.ind~party.comp+ executive.comp + military, data=a.out, model="ls") 
m71<-zelig(preffalse.unadj.ind~party.comp+ executive.comp  +military + repression + duration, data=a.out, model="ls") 
m72<-zelig(preffalse.unadj.ind~party.comp+ executive.comp + military + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m73<-zelig(preffalse.unadj.ind~repression, data=a.out, model="ls") 
m74<-zelig(preffalse.unadj.ind~repression+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m75<-zelig(preffalse.unadj.ind~repression + party.comp+ executive.comp  +military +duration, data=a.out, model="ls") 
m76<-zelig(preffalse.ind~repression+ party.comp+ executive.comp + military  + duration + oil + gdppc + edu, data=a.out, model="ls") 
m77<-zelig(preffalse.unadj.ind~duration, data=a.out, model="ls") 
m78<-zelig(preffalse.unadj.ind~duration+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m79<-zelig(preffalse.unadj.ind~duration+ party.comp+ executive.comp + military + repression, data=a.out, model="ls") 
m80<-zelig(preffalse.unadj.ind~duration +party.comp+ executive.comp + military + repression + oil + gdppc + edu, data=a.out, model="ls") 
m81<-zelig(preffalse.unadj.ind.min~executive.comp, data=a.out, model="ls") 
m82<-zelig(preffalse.unadj.ind.min~executive.comp + military + party.comp, data=a.out, model="ls") 
m83<-zelig(preffalse.unadj.ind.min~executive.comp + military + party.comp + repression + duration, data=a.out, model="ls") 
m84<-zelig(preffalse.unadj.ind.min~executive.comp + military + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m85<-zelig(preffalse.unadj.ind.min~military, data=a.out, model="ls") 
m86<-zelig(preffalse.unadj.ind.min~military+ executive.comp + party.comp, data=a.out, model="ls") 
m87<-zelig(preffalse.unadj.ind.min~military+ executive.comp  + party.comp + repression + duration, data=a.out, model="ls") 
m88<-zelig(preffalse.unadj.ind.min~military+ executive.comp + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m89<-zelig(preffalse.unadj.ind.min~party.comp, data=a.out, model="ls") 
m90<-zelig(preffalse.unadj.ind.min~party.comp+ executive.comp + military, data=a.out, model="ls") 
m91<-zelig(preffalse.unadj.ind.min~party.comp+ executive.comp  +military + repression + duration, data=a.out, model="ls") 
m92<-zelig(preffalse.unadj.ind.min~party.comp+ executive.comp + military + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m93<-zelig(preffalse.unadj.ind.min~repression, data=a.out, model="ls") 
m94<-zelig(preffalse.unadj.ind.min~repression+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m95<-zelig(preffalse.unadj.ind.min~repression + party.comp+ executive.comp  +military +duration, data=a.out, model="ls") 
m96<-zelig(preffalse.unadj.ind.min~repression+ party.comp+ executive.comp + military  + duration + oil + gdppc + edu, data=a.out, model="ls") 
m97<-zelig(preffalse.unadj.ind.min~duration, data=a.out, model="ls") 
m98<-zelig(preffalse.unadj.ind.min~duration+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m99<-zelig(preffalse.unadj.ind.min~duration+ party.comp+ executive.comp  +military + repression, data=a.out, model="ls") 
m100<-zelig(preffalse.unadj.ind.min~duration +party.comp+ executive.comp + military + repression + oil + gdppc + edu, data=a.out, model="ls") 
m101<-zelig(preffalse.unadj.ind.mid~executive.comp, data=a.out, model="ls") 
m102<-zelig(preffalse.unadj.ind.mid~executive.comp + military + party.comp, data=a.out, model="ls") 
m103<-zelig(preffalse.unadj.ind.mid~executive.comp + military + party.comp + repression + duration, data=a.out, model="ls") 
m104<-zelig(preffalse.unadj.ind.mid~executive.comp + military + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m105<-zelig(preffalse.unadj.ind.mid~military, data=a.out, model="ls") 
m106<-zelig(preffalse.unadj.ind.mid~military+ executive.comp + party.comp, data=a.out, model="ls") 
m107<-zelig(preffalse.unadj.ind.mid~military+ executive.comp  + party.comp + repression + duration, data=a.out, model="ls") 
m108<-zelig(preffalse.unadj.ind.mid~military+ executive.comp + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m109<-zelig(preffalse.unadj.ind.mid~party.comp, data=a.out, model="ls") 
m110<-zelig(preffalse.unadj.ind.mid~party.comp+ executive.comp + military, data=a.out, model="ls") 
m111<-zelig(preffalse.unadj.ind.mid~party.comp+ executive.comp  +military + repression + duration, data=a.out, model="ls") 
m112<-zelig(preffalse.unadj.ind.mid~party.comp+ executive.comp + military + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m113<-zelig(preffalse.unadj.ind.mid~repression, data=a.out, model="ls") 
m114<-zelig(preffalse.unadj.ind.mid~repression+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m115<-zelig(preffalse.unadj.ind.mid~repression + party.comp+ executive.comp  +military +duration, data=a.out, model="ls") 
m116<-zelig(preffalse.unadj.ind.mid~repression+ party.comp+ executive.comp + military  + duration + oil + gdppc + edu, data=a.out, model="ls") 
m117<-zelig(preffalse.unadj.ind.mid~duration, data=a.out, model="ls") 
m118<-zelig(preffalse.unadj.ind.mid~duration+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m119<-zelig(preffalse.unadj.ind.mid~duration+ party.comp+ executive.comp  +military + repression, data=a.out, model="ls") 
m120<-zelig(preffalse.unadj.ind.mid~duration +party.comp+ executive.comp + military + repression + oil + gdppc + edu, data=a.out, model="ls")
m121<-zelig(preffalse.ind.20~executive.comp, data=a.out, model="ls") 
m122<-zelig(preffalse.ind.20~executive.comp + military + party.comp, data=a.out, model="ls") 
m123<-zelig(preffalse.ind.20~executive.comp + military + party.comp + repression + duration, data=a.out, model="ls") 
m124<-zelig(preffalse.ind.20~executive.comp + military + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m125<-zelig(preffalse.ind.20~military, data=a.out, model="ls") 
m126<-zelig(preffalse.ind.20~military+ executive.comp + party.comp, data=a.out, model="ls") 
m127<-zelig(preffalse.ind.20~military+ executive.comp  + party.comp + repression + duration, data=a.out, model="ls") 
m128<-zelig(preffalse.ind.20~military+ executive.comp + party.comp + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m129<-zelig(preffalse.ind.20~party.comp, data=a.out, model="ls") 
m130<-zelig(preffalse.ind.20~party.comp+ executive.comp + military, data=a.out, model="ls") 
m131<-zelig(preffalse.ind.20~party.comp+ executive.comp  +military + repression + duration, data=a.out, model="ls") 
m132<-zelig(preffalse.ind.20~party.comp+ executive.comp + military + repression + duration + oil + gdppc + edu, data=a.out, model="ls") 
m133<-zelig(preffalse.ind.20~repression, data=a.out, model="ls") 
m134<-zelig(preffalse.ind.20~repression+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m135<-zelig(preffalse.ind.20~repression + party.comp+ executive.comp  +military +duration, data=a.out, model="ls") 
m136<-zelig(preffalse.ind.20~repression+ party.comp+ executive.comp + military  + duration + oil + gdppc + edu, data=a.out, model="ls") 
m137<-zelig(preffalse.ind.20~duration, data=a.out, model="ls") 
m138<-zelig(preffalse.ind.20~duration+ party.comp+ executive.comp + military, data=a.out, model="ls") 
m139<-zelig(preffalse.ind.20~duration+ party.comp+ executive.comp  +military + repression, data=a.out, model="ls") 
m140<-zelig(preffalse.ind.20~duration +party.comp+ executive.comp + military + repression + oil + gdppc + edu, data=a.out, model="ls") 

ispf.models<-read.csv("ispf.models.csv",stringsAsFactors=FALSE) 
for (j in 1:140) {
  try(eval(parse(text=paste("ispf.models$coeff[",j,"]<-unlist(combine_coef_se(m",j,"))[2,1]",sep=""))))
  try(eval(parse(text=paste("ispf.models$se[",j,"]<-unlist(combine_coef_se(m",j,"))[2,2]",sep=""))))
  try(eval(parse(text=paste("ispf.models$pvalue[",j,"]<-unlist(combine_coef_se(m",j,"))[2,4]",sep=""))))
}  
write.csv(ispf.models, "ispf.models.results.csv")

ispf.models$l95ci<-ispf.models$coeff-1.96*ispf.models$se
ispf.models$u95ci<-ispf.models$coeff+1.96*ispf.models$se

ispf.models<-subset(ispf.models, ispf.models$include==1)

ispf.models.executive.comp<-subset(ispf.models, ispf.models$iv=="executive.comp")
ispf.models.military<-subset(ispf.models, ispf.models$iv=="military")
ispf.models.party.comp<-subset(ispf.models, ispf.models$iv=="party.comp")
ispf.models.repression<-subset(ispf.models, ispf.models$iv=="repression")
ispf.models.duration<-subset(ispf.models, ispf.models$iv=="duration")

pdf('fig-robustness-adj.pdf', width=9.5, height=10) 
p1<-ggplot(ispf.models.executive.comp, aes(x=pvalue, y=coeff, shape=type,color=dv)) + xlab("p-value") + ylab("Coefficient Estimate") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.85) + ggtitle("Executive Competition (exec.comp)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color="grey50") + geom_vline(xintercept =.1,linetype =3,color="grey50") + geom_vline(xintercept =.01,linetype =3,color="grey50") + xlim(0, 1) + ylim(-.06, .06) + scale_color_manual(values=c("red4","red1","rosybrown","dodgerblue4"), labels = c("cens.ind.3q.3q", "cens.ind.3q.20q", "cens.ind.1q.1q", "cens.ind.3q.0q"), name="Dependent Variable") + scale_shape_manual(values=c(0,1,2,5), breaks = c("bivariate", "institutions","repression + duration","full"), name="Covariate Set") + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position="right") + geom_hline(yintercept =.0,linetype =2,color="grey50") + theme(legend.key.height = unit(.45, "cm"))
p2<-ggplot(ispf.models.military, aes(x=pvalue, y=coeff, shape=type,color=dv)) + xlab("p-value") + ylab("Coefficient Estimate") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.85) + ggtitle("Military Control (military)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color="grey50") + geom_vline(xintercept =.1,linetype =3,color="grey50") + geom_vline(xintercept =.01,linetype =3,color="grey50") + xlim(0, 1) + ylim(-.06, .06) + scale_color_manual(values=c("red4","red1","rosybrown","dodgerblue4"), labels = c("cens.ind.3q.3q", "cens.ind.3q.20q", "cens.ind.1q.1q", "cens.ind.3q.0q"), name="Dependent Variable") + scale_shape_manual(values=c(0,1,2,5), breaks = c("bivariate", "institutions","repression + duration","full"), name="Covariate Set") + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position="right") + geom_hline(yintercept =.0,linetype =2,color="grey50") + theme(legend.key.height = unit(.45, "cm"))
p3<-ggplot(ispf.models.party.comp, aes(x=pvalue, y=coeff, shape=type,color=dv)) + xlab("p-value") + ylab("Coefficient Estimate") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.85) + ggtitle("Party Competition (party.comp)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color="grey50") + geom_vline(xintercept =.1,linetype =3,color="grey50") + geom_vline(xintercept =.01,linetype =3,color="grey50") + xlim(0, 1) + ylim(-.06, .06) + scale_color_manual(values=c("red4","red1","rosybrown","dodgerblue4"), labels = c("cens.ind.3q.3q", "cens.ind.3q.20q", "cens.ind.1q.1q", "cens.ind.3q.0q"), name="Dependent Variable") + scale_shape_manual(values=c(0,1,2,5), breaks = c("bivariate", "institutions","repression + duration","full"), name="Covariate Set") + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position="right") + geom_hline(yintercept =.0,linetype =2,color="grey50") + theme(legend.key.height = unit(.45, "cm"))
p4<-ggplot(ispf.models.repression, aes(x=pvalue, y=coeff, shape=type,color=dv)) + xlab("p-value") + ylab("Coefficient Estimate") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.85) + ggtitle("Repression (repression)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color="grey50") + geom_vline(xintercept =.1,linetype =3,color="grey50") + geom_vline(xintercept =.01,linetype =3,color="grey50") + xlim(0, 1) + ylim(-.06, .06) + scale_color_manual(values=c("red4","red1","rosybrown","dodgerblue4"), labels = c("cens.ind.3q.3q", "cens.ind.3q.20q", "cens.ind.1q.1q", "cens.ind.3q.0q"), name="Dependent Variable") + scale_shape_manual(values=c(0,1,2,5), breaks = c("bivariate", "institutions","repression + duration","full"), name="Covariate Set") + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position="right") + geom_hline(yintercept =.0,linetype =2,color="grey50") + theme(legend.key.height = unit(.45, "cm"))
p5<-ggplot(ispf.models.duration, aes(x=pvalue, y=coeff, shape=type,color=dv)) + xlab("p-value") + ylab("Coefficient Estimate") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.85) + ggtitle("Regime Duration (duration)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) +geom_vline(xintercept =.05,linetype =3,color="grey50") + geom_vline(xintercept =.1,linetype =3,color="grey50") + geom_vline(xintercept =.01,linetype =3,color="grey50") + xlim(0, 1) + ylim(-.06, .06) + scale_color_manual(values=c("red4","red1","rosybrown","dodgerblue4"), labels = c("cens.ind.3q.3q", "cens.ind.3q.20q", "cens.ind.1q.1q", "cens.ind.3q.0q"), name="Dependent Variable") + scale_shape_manual(values=c(0,1,2,5), breaks = c("bivariate", "institutions","repression + duration","full"), name="Covariate Set") + theme(plot.title = element_text(size = 12),legend.title=element_text(size=8), legend.text=element_text(size=6)) + theme(legend.position="right") + geom_hline(yintercept =.0,linetype =2,color="grey50") + theme(legend.key.height = unit(.45, "cm"))
arrange_ggplot2(p1,p3,p5,p2,p4, nrow=3, ncol=2)
dev.off()


###INDEX PERFORMANCE###

cor.test(data.impute$E069_11.dkna.rate, data.impute$E124.dkna.rate, method = "pearson", conf.level = 0.95)
cor.test(data.impute$E069_11.dkna.rate, data.impute$E236.dkna.rate, method = "pearson", conf.level = 0.95)
cor.test(data.impute$E124.dkna.rate, data.impute$E236.dkna.rate, method = "pearson", conf.level = 0.95)

pdf('fig-correlations-sens.pdf', width=9.5, height=3.25)
p1<-ggplot(data.impute, aes(x=E069_11.dkna.rate, y=E124.dkna.rate)) + xlab("conf.gov (item nonreponse rate)") + ylab("humanrights (item nonresponse rate)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.543")+xlim(c(0,.25))+ylim(c(0,.25))
p2<-ggplot(data.impute, aes(x=E069_11.dkna.rate, y=E236.dkna.rate)) + xlab("conf.gov (item nonreponse rate)") + ylab("democraticness (item nonresponse rate)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.672")+xlim(c(0,.25))+ylim(c(0,.25))
p3<-ggplot(data.impute, aes(x=E124.dkna.rate, y=E236.dkna.rate)) + xlab("humanrights (item nonresponse rate)") + ylab("democraticness (item nonresponse rate)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.939")+xlim(c(0,.25))+ylim(c(0,.25))
arrange_ggplot2(p1,p2,p3, nrow=1, ncol=3)
dev.off()

cor.test(data.impute$E069_10.dkna.rate, data.impute$A170.dkna.rate, method = "pearson", conf.level = 0.95)
cor.test(data.impute$E069_10.dkna.rate, data.impute$A165.dkna.rate, method = "pearson", conf.level = 0.95)
cor.test(data.impute$A170.dkna.rate, data.impute$A165.dkna.rate, method = "pearson", conf.level = 0.95)

pdf('fig-correlations-nonsens.pdf', width=9.5, height=3.25)
p1<-ggplot(data.impute, aes(x=E069_10.dkna.rate, y=A170.dkna.rate)) + xlab("conf.tv (item nonreponse rate)") + ylab("life.sat (item nonresponse rate)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.486")+xlim(c(0,.25))+ylim(c(0,.25))
p2<-ggplot(data.impute, aes(x=E069_10.dkna.rate, y=A165.dkna.rate)) + xlab("conf.tv (item nonreponse rate)") + ylab("gen.trust (item nonresponse rate)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.319")+xlim(c(0,.25))+ylim(c(0,.25))
p3<-ggplot(data.impute, aes(x=A170.dkna.rate, y=A165.dkna.rate)) + xlab("life.sat (item nonresponse rate)") + ylab("gen.trust (item nonresponse rate)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.283")+xlim(c(0,.25))+ylim(c(0,.25))
arrange_ggplot2(p1,p2,p3, nrow=1, ncol=3)
dev.off()

cor.test(outdata$preffalse.ind, outdata$preffalse.ind.20, method = "pearson", conf.level = 0.95)
cor.test(outdata$preffalse.ind, outdata$preffalse.ind.min, method = "pearson", conf.level = 0.95)
cor.test(outdata$preffalse.ind, outdata$preffalse.unadj.ind, method = "pearson", conf.level = 0.95)

pdf('fig-correlations-preffalse.pdf', width=9.5, height=3.25)
p1<-ggplot(outdata, aes(x=preffalse.ind, y=preffalse.ind.20)) + xlab("cens.ind.3q.3q") + ylab("cens.ind.3q.20q") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.762")+xlim(c(-.05,.25))+ylim(c(-0.05,.25))
p2<-ggplot(outdata, aes(x=preffalse.ind, y=preffalse.ind.min)) + xlab("cens.ind.3q.3q") + ylab("cens.ind.1q.1q") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.183")+xlim(c(-0.05,.25))+ylim(c(-0.05,.25))
p3<-ggplot(outdata, aes(x=preffalse.ind, y=preffalse.unadj.ind)) + xlab("cens.ind.3q.3q") + ylab("cens.ind.3q.0q") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_smooth(method='lm', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .22, y = .25, label = "r = 0.825")+xlim(c(-0.05,.25))+ylim(c(-0.05,.25))
arrange_ggplot2(p1,p2,p3, nrow=1, ncol=3)
dev.off() 

cor.test(outdata$nonsens.dkna.rate, outdata$sens.dkna.rate, method = "pearson", conf.level = 0.95)

pdf('fig-correlations-sensnonsens.pdf', width=5, height=5)
ggplot(outdata, aes(x=nonsens.dkna.rate, y=sens.dkna.rate)) + xlab("Sample Item Nonresponse Rate (Nonsensitive Questions)") + ylab("Sample Item Nonresponse Rate (Sensitive Questions)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + theme_bw() + geom_point(alpha=.5)  + theme(plot.title = element_text(hjust = 0.5)) + geom_abline(intercept = 0, slope = 1) +xlim(c(0,.28))+ylim(c(0,.28))  + geom_smooth(method='loess', se=FALSE, color="red", linetype="dashed",size=.5) +  annotate("text", x = .27, y = .28, label = "r = 0.776")
dev.off()

#####CHINA ANALYSIS#####

#WVS - Wave 5
wvs.china<-subset(wvs, wvs$S003=="China")
wvs.china.w5<-subset(wvs.china, wvs.china$S025=="China (2007)")

wvs.china.w5$E069_11
wvs.china.w5$conf.gov<-NA
wvs.china.w5$conf.gov[wvs.china.w5$E069_11=="A great deal"]<-4
wvs.china.w5$conf.gov[wvs.china.w5$E069_11=="Quite a lot"]<-3
wvs.china.w5$conf.gov[wvs.china.w5$E069_11=="Not very much"]<-2
wvs.china.w5$conf.gov[wvs.china.w5$E069_11=="None at all"]<-1

wvs.china.w5$birthyear<-NA
wvs.china.w5$birthyear[as.numeric(paste(wvs.china.w5$X002))<2000]<-"1980+"
wvs.china.w5$birthyear[as.numeric(paste(wvs.china.w5$X002))<1980]<-"1970s"
wvs.china.w5$birthyear[as.numeric(paste(wvs.china.w5$X002))<1970]<-"1960s"
wvs.china.w5$birthyear[as.numeric(paste(wvs.china.w5$X002))<1960]<-"1950s"
wvs.china.w5$birthyear[as.numeric(paste(wvs.china.w5$X002))<1950]<-"<1950"

wvs.china.w5$institution<-NA
wvs.china.w5$institution[wvs.china.w5$X052=="Public institution"]<-"Public"
wvs.china.w5$institution[wvs.china.w5$X052=="Self-employed"]<-"Private"
wvs.china.w5$institution[wvs.china.w5$X052=="Private business"]<-"Private"
wvs.china.w5$institution[wvs.china.w5$X052=="Private non-profit organization"]<-"Private"

wvs.china.w5$education<-NA
wvs.china.w5$education[wvs.china.w5$X025=="Complete secondary school: technical/vocational type/Seconda"]<-"Secondary"
wvs.china.w5$education[wvs.china.w5$X025=="Complete secondary: university-preparatory type/Full seconda"]<-"Secondary"
wvs.china.w5$education[wvs.china.w5$X025=="Completed (compulsory) elementary education"]<-"Elementary"
wvs.china.w5$education[wvs.china.w5$X025=="University with degree/Higher education - upper-level tertia"]<-"University"
wvs.china.w5$education[wvs.china.w5$X025=="Not applicable; No formal education"]<-"No education"

wvs.china.w5$class<-NA
wvs.china.w5$class[wvs.china.w5$X045=="Lower class"]<-"1 - Lower class"
wvs.china.w5$class[wvs.china.w5$X045=="Working class"]<-"2 - Lower middle class"
wvs.china.w5$class[wvs.china.w5$X045=="Lower middle class"]<-"3 - Middle class"
wvs.china.w5$class[wvs.china.w5$X045=="Upper middle class"]<-"4 - Upper middle class"
wvs.china.w5$class[wvs.china.w5$X045=="Upper class"]<-"5 - Upper class"

data.china.w5.age<-aggregate(wvs.china.w5$conf.gov.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)
data.china.w5.age[,3]<-aggregate(wvs.china.w5$conf.gov.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,4]<-aggregate(wvs.china.w5$resp.hr.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,5]<-aggregate(wvs.china.w5$resp.hr.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,6]<-aggregate(wvs.china.w5$democraticness.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,7]<-aggregate(wvs.china.w5$democraticness.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,8]<-aggregate(wvs.china.w5$conf.tv.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,9]<-aggregate(wvs.china.w5$conf.tv.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,10]<-aggregate(wvs.china.w5$gentrust.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,11]<-aggregate(wvs.china.w5$gentrust.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,12]<-aggregate(wvs.china.w5$lifesat.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,13]<-aggregate(wvs.china.w5$lifesat.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,14]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$birthyear), FUN="mean", na.rm=TRUE)[2]
data.china.w5.age[,15]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$birthyear), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w5.age)<-c("birthyear","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w5.age$conf.gov.dk.rate<-data.china.w5.age$conf.gov.dk/data.china.w5.age$conf.gov.ob
data.china.w5.age$resp.hr.dk.rate<-data.china.w5.age$resp.hr.dk/data.china.w5.age$resp.hr.ob
data.china.w5.age$democraticness.dk.rate<-data.china.w5.age$democraticness.dk/data.china.w5.age$democraticness.ob
data.china.w5.age$conf.tv.dk.rate<-data.china.w5.age$conf.tv.dk/data.china.w5.age$conf.tv.ob
data.china.w5.age$gentrust.dk.rate<-data.china.w5.age$gentrust.dk/data.china.w5.age$gentrust.ob
data.china.w5.age$lifesat.dk.rate<-data.china.w5.age$lifesat.dk/data.china.w5.age$lifesat.ob
data.china.w5.age$sens.ob<-data.china.w5.age$conf.gov.ob+data.china.w5.age$resp.hr.ob+data.china.w5.age$democraticness.ob
data.china.w5.age$sens.ob[data.china.w5.age$conf.gov.ob==0]<-NA
data.china.w5.age$sens.ob[data.china.w5.age$resp.hr.ob==0]<-NA
data.china.w5.age$sens.ob[data.china.w5.age$democraticness.ob==0]<-NA
data.china.w5.age$sens.dk<-data.china.w5.age$conf.gov.dk+data.china.w5.age$resp.hr.dk+data.china.w5.age$democraticness.dk
data.china.w5.age$sens.dk[data.china.w5.age$conf.gov.ob==0]<-NA
data.china.w5.age$sens.dk[data.china.w5.age$resp.hr.ob==0]<-NA
data.china.w5.age$sens.dk[data.china.w5.age$democraticness.ob==0]<-NA
data.china.w5.age$sens.dk.rate<-data.china.w5.age$sens.dk/data.china.w5.age$sens.ob
data.china.w5.age$sens.dk.se<-sqrt((data.china.w5.age$sens.dk.rate*(1-data.china.w5.age$sens.dk.rate))/data.china.w5.age$sens.ob)
data.china.w5.age$sens.dk.u95ci<-data.china.w5.age$sens.dk.rate+1.96*data.china.w5.age$sens.dk.se
data.china.w5.age$sens.dk.l95ci<-data.china.w5.age$sens.dk.rate-1.96*data.china.w5.age$sens.dk.se
data.china.w5.age$nonsens.ob<-data.china.w5.age$conf.tv.ob+data.china.w5.age$lifesat.ob+data.china.w5.age$gentrust.ob
data.china.w5.age$nonsens.ob[data.china.w5.age$conf.tv.ob==0]<-NA
data.china.w5.age$nonsens.ob[data.china.w5.age$lifesat.ob==0]<-NA
data.china.w5.age$nonsens.ob[data.china.w5.age$gentrust.ob==0]<-NA
data.china.w5.age$nonsens.dk<-data.china.w5.age$conf.tv.dk+data.china.w5.age$lifesat.dk+data.china.w5.age$gentrust.dk
data.china.w5.age$nonsens.dk[data.china.w5.age$conf.tv.ob==0]<-NA
data.china.w5.age$nonsens.dk[data.china.w5.age$lifesat.ob==0]<-NA
data.china.w5.age$nonsens.dk[data.china.w5.age$gentrust.ob==0]<-NA
data.china.w5.age$nonsens.dk.rate<-data.china.w5.age$nonsens.dk/data.china.w5.age$nonsens.ob
data.china.w5.age$nonsens.dk.se<-sqrt((data.china.w5.age$nonsens.dk.rate*(1-data.china.w5.age$nonsens.dk.rate))/data.china.w5.age$nonsens.ob)
data.china.w5.age$nonsens.dk.u95ci<-data.china.w5.age$nonsens.dk.rate+1.96*data.china.w5.age$nonsens.dk.se
data.china.w5.age$nonsens.dk.l95ci<-data.china.w5.age$nonsens.dk.rate-1.96*data.china.w5.age$nonsens.dk.se
data.china.w5.age$preffalse.ind<-data.china.w5.age$sens.dk.rate-data.china.w5.age$nonsens.dk.rate
data.china.w5.age$preffalse.ind.se<-sqrt((data.china.w5.age$nonsens.dk.se^2)+(data.china.w5.age$sens.dk.se^2))
data.china.w5.age$preffalse.ind.u95ci<-data.china.w5.age$preffalse.ind+1.96*data.china.w5.age$preffalse.ind.se
data.china.w5.age$preffalse.ind.l95ci<-data.china.w5.age$preffalse.ind-1.96*data.china.w5.age$preffalse.ind.se
data.china.w5.age$conf.gov.se<-data.china.w5.age$conf.gov.sd/sqrt(data.china.w5.age$conf.gov.ob)
data.china.w5.age$conf.gov.u95ci<-data.china.w5.age$conf.gov.mean+1.96*data.china.w5.age$conf.gov.se
data.china.w5.age$conf.gov.l95ci<-data.china.w5.age$conf.gov.mean-1.96*data.china.w5.age$conf.gov.se
data.china.w5.age$variable<-"Birthyear"
colnames(data.china.w5.age)[1] <- "level"

data.china.w5.gender<-aggregate(wvs.china.w5$conf.gov.ob, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)
data.china.w5.gender[,3]<-aggregate(wvs.china.w5$conf.gov.dk, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,4]<-aggregate(wvs.china.w5$resp.hr.ob, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,5]<-aggregate(wvs.china.w5$resp.hr.dk, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,6]<-aggregate(wvs.china.w5$democraticness.ob, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,7]<-aggregate(wvs.china.w5$democraticness.dk, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,8]<-aggregate(wvs.china.w5$conf.tv.ob, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,9]<-aggregate(wvs.china.w5$conf.tv.dk, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,10]<-aggregate(wvs.china.w5$gentrust.ob, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,11]<-aggregate(wvs.china.w5$gentrust.dk, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,12]<-aggregate(wvs.china.w5$lifesat.ob, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,13]<-aggregate(wvs.china.w5$lifesat.dk, by=list(wvs.china.w5$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w5.gender[,14]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$X001), FUN="mean", na.rm=TRUE)[2]
data.china.w5.gender[,15]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$X001), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w5.gender)<-c("gender","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w5.gender$conf.gov.dk.rate<-data.china.w5.gender$conf.gov.dk/data.china.w5.gender$conf.gov.ob
data.china.w5.gender$resp.hr.dk.rate<-data.china.w5.gender$resp.hr.dk/data.china.w5.gender$resp.hr.ob
data.china.w5.gender$democraticness.dk.rate<-data.china.w5.gender$democraticness.dk/data.china.w5.gender$democraticness.ob
data.china.w5.gender$conf.tv.dk.rate<-data.china.w5.gender$conf.tv.dk/data.china.w5.gender$conf.tv.ob
data.china.w5.gender$gentrust.dk.rate<-data.china.w5.gender$gentrust.dk/data.china.w5.gender$gentrust.ob
data.china.w5.gender$lifesat.dk.rate<-data.china.w5.gender$lifesat.dk/data.china.w5.gender$lifesat.ob
data.china.w5.gender$sens.ob<-data.china.w5.gender$conf.gov.ob+data.china.w5.gender$resp.hr.ob+data.china.w5.gender$democraticness.ob
data.china.w5.gender$sens.ob[data.china.w5.gender$conf.gov.ob==0]<-NA
data.china.w5.gender$sens.ob[data.china.w5.gender$resp.hr.ob==0]<-NA
data.china.w5.gender$sens.ob[data.china.w5.gender$democraticness.ob==0]<-NA
data.china.w5.gender$sens.dk<-data.china.w5.gender$conf.gov.dk+data.china.w5.gender$resp.hr.dk+data.china.w5.gender$democraticness.dk
data.china.w5.gender$sens.dk[data.china.w5.gender$conf.gov.ob==0]<-NA
data.china.w5.gender$sens.dk[data.china.w5.gender$resp.hr.ob==0]<-NA
data.china.w5.gender$sens.dk[data.china.w5.gender$democraticness.ob==0]<-NA
data.china.w5.gender$sens.dk.rate<-data.china.w5.gender$sens.dk/data.china.w5.gender$sens.ob
data.china.w5.gender$sens.dk.se<-sqrt((data.china.w5.gender$sens.dk.rate*(1-data.china.w5.gender$sens.dk.rate))/data.china.w5.gender$sens.ob)
data.china.w5.gender$sens.dk.u95ci<-data.china.w5.gender$sens.dk.rate+1.96*data.china.w5.gender$sens.dk.se
data.china.w5.gender$sens.dk.l95ci<-data.china.w5.gender$sens.dk.rate-1.96*data.china.w5.gender$sens.dk.se
data.china.w5.gender$nonsens.ob<-data.china.w5.gender$conf.tv.ob+data.china.w5.gender$lifesat.ob+data.china.w5.gender$gentrust.ob
data.china.w5.gender$nonsens.ob[data.china.w5.gender$conf.tv.ob==0]<-NA
data.china.w5.gender$nonsens.ob[data.china.w5.gender$lifesat.ob==0]<-NA
data.china.w5.gender$nonsens.ob[data.china.w5.gender$gentrust.ob==0]<-NA
data.china.w5.gender$nonsens.dk<-data.china.w5.gender$conf.tv.dk+data.china.w5.gender$lifesat.dk+data.china.w5.gender$gentrust.dk
data.china.w5.gender$nonsens.dk[data.china.w5.gender$conf.tv.ob==0]<-NA
data.china.w5.gender$nonsens.dk[data.china.w5.gender$lifesat.ob==0]<-NA
data.china.w5.gender$nonsens.dk[data.china.w5.gender$gentrust.ob==0]<-NA
data.china.w5.gender$nonsens.dk.rate<-data.china.w5.gender$nonsens.dk/data.china.w5.gender$nonsens.ob
data.china.w5.gender$nonsens.dk.se<-sqrt((data.china.w5.gender$nonsens.dk.rate*(1-data.china.w5.gender$nonsens.dk.rate))/data.china.w5.gender$nonsens.ob)
data.china.w5.gender$nonsens.dk.u95ci<-data.china.w5.gender$nonsens.dk.rate+1.96*data.china.w5.gender$nonsens.dk.se
data.china.w5.gender$nonsens.dk.l95ci<-data.china.w5.gender$nonsens.dk.rate-1.96*data.china.w5.gender$nonsens.dk.se
data.china.w5.gender$preffalse.ind<-data.china.w5.gender$sens.dk.rate-data.china.w5.gender$nonsens.dk.rate
data.china.w5.gender$preffalse.ind.se<-sqrt((data.china.w5.gender$nonsens.dk.se^2)+(data.china.w5.gender$sens.dk.se^2))
data.china.w5.gender$preffalse.ind.u95ci<-data.china.w5.gender$preffalse.ind+1.96*data.china.w5.gender$preffalse.ind.se
data.china.w5.gender$preffalse.ind.l95ci<-data.china.w5.gender$preffalse.ind-1.96*data.china.w5.gender$preffalse.ind.se
data.china.w5.gender$conf.gov.se<-data.china.w5.gender$conf.gov.sd/sqrt(data.china.w5.gender$conf.gov.ob)
data.china.w5.gender$conf.gov.u95ci<-data.china.w5.gender$conf.gov.mean+1.96*data.china.w5.gender$conf.gov.se
data.china.w5.gender$conf.gov.l95ci<-data.china.w5.gender$conf.gov.mean-1.96*data.china.w5.gender$conf.gov.se
ggplot(data.china.w5.gender, aes(x=gender, y=preffalse.ind)) + xlab("Gender") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5)  +  geom_segment(aes(x = gender, y = preffalse.ind.l95ci, xend = gender, yend = preffalse.ind.u95ci), color="grey60", alpha=.5, data = data.china.w5.gender) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + geom_hline(aes(yintercept=0.00), colour="grey40", linetype="dotted") 
data.china.w5.gender$variable<-"Gender"
colnames(data.china.w5.gender)[1] <- "level"

data.china.w5.education<-aggregate(wvs.china.w5$conf.gov.ob, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)
data.china.w5.education[,3]<-aggregate(wvs.china.w5$conf.gov.dk, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,4]<-aggregate(wvs.china.w5$resp.hr.ob, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,5]<-aggregate(wvs.china.w5$resp.hr.dk, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,6]<-aggregate(wvs.china.w5$democraticness.ob, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,7]<-aggregate(wvs.china.w5$democraticness.dk, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,8]<-aggregate(wvs.china.w5$conf.tv.ob, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,9]<-aggregate(wvs.china.w5$conf.tv.dk, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,10]<-aggregate(wvs.china.w5$gentrust.ob, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,11]<-aggregate(wvs.china.w5$gentrust.dk, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,12]<-aggregate(wvs.china.w5$lifesat.ob, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,13]<-aggregate(wvs.china.w5$lifesat.dk, by=list(wvs.china.w5$education), FUN="sum", na.rm=TRUE)[2]
data.china.w5.education[,14]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$education), FUN="mean", na.rm=TRUE)[2]
data.china.w5.education[,15]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$education), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w5.education)<-c("education","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w5.education$conf.gov.dk.rate<-data.china.w5.education$conf.gov.dk/data.china.w5.education$conf.gov.ob
data.china.w5.education$resp.hr.dk.rate<-data.china.w5.education$resp.hr.dk/data.china.w5.education$resp.hr.ob
data.china.w5.education$democraticness.dk.rate<-data.china.w5.education$democraticness.dk/data.china.w5.education$democraticness.ob
data.china.w5.education$conf.tv.dk.rate<-data.china.w5.education$conf.tv.dk/data.china.w5.education$conf.tv.ob
data.china.w5.education$gentrust.dk.rate<-data.china.w5.education$gentrust.dk/data.china.w5.education$gentrust.ob
data.china.w5.education$lifesat.dk.rate<-data.china.w5.education$lifesat.dk/data.china.w5.education$lifesat.ob
data.china.w5.education$sens.ob<-data.china.w5.education$conf.gov.ob+data.china.w5.education$resp.hr.ob+data.china.w5.education$democraticness.ob
data.china.w5.education$sens.ob[data.china.w5.education$conf.gov.ob==0]<-NA
data.china.w5.education$sens.ob[data.china.w5.education$resp.hr.ob==0]<-NA
data.china.w5.education$sens.ob[data.china.w5.education$democraticness.ob==0]<-NA
data.china.w5.education$sens.dk<-data.china.w5.education$conf.gov.dk+data.china.w5.education$resp.hr.dk+data.china.w5.education$democraticness.dk
data.china.w5.education$sens.dk[data.china.w5.education$conf.gov.ob==0]<-NA
data.china.w5.education$sens.dk[data.china.w5.education$resp.hr.ob==0]<-NA
data.china.w5.education$sens.dk[data.china.w5.education$democraticness.ob==0]<-NA
data.china.w5.education$sens.dk.rate<-data.china.w5.education$sens.dk/data.china.w5.education$sens.ob
data.china.w5.education$sens.dk.se<-sqrt((data.china.w5.education$sens.dk.rate*(1-data.china.w5.education$sens.dk.rate))/data.china.w5.education$sens.ob)
data.china.w5.education$sens.dk.u95ci<-data.china.w5.education$sens.dk.rate+1.96*data.china.w5.education$sens.dk.se
data.china.w5.education$sens.dk.l95ci<-data.china.w5.education$sens.dk.rate-1.96*data.china.w5.education$sens.dk.se
data.china.w5.education$nonsens.ob<-data.china.w5.education$conf.tv.ob+data.china.w5.education$lifesat.ob+data.china.w5.education$gentrust.ob
data.china.w5.education$nonsens.ob[data.china.w5.education$conf.tv.ob==0]<-NA
data.china.w5.education$nonsens.ob[data.china.w5.education$lifesat.ob==0]<-NA
data.china.w5.education$nonsens.ob[data.china.w5.education$gentrust.ob==0]<-NA
data.china.w5.education$nonsens.dk<-data.china.w5.education$conf.tv.dk+data.china.w5.education$lifesat.dk+data.china.w5.education$gentrust.dk
data.china.w5.education$nonsens.dk[data.china.w5.education$conf.tv.ob==0]<-NA
data.china.w5.education$nonsens.dk[data.china.w5.education$lifesat.ob==0]<-NA
data.china.w5.education$nonsens.dk[data.china.w5.education$gentrust.ob==0]<-NA
data.china.w5.education$nonsens.dk.rate<-data.china.w5.education$nonsens.dk/data.china.w5.education$nonsens.ob
data.china.w5.education$nonsens.dk.se<-sqrt((data.china.w5.education$nonsens.dk.rate*(1-data.china.w5.education$nonsens.dk.rate))/data.china.w5.education$nonsens.ob)
data.china.w5.education$nonsens.dk.u95ci<-data.china.w5.education$nonsens.dk.rate+1.96*data.china.w5.education$nonsens.dk.se
data.china.w5.education$nonsens.dk.l95ci<-data.china.w5.education$nonsens.dk.rate-1.96*data.china.w5.education$nonsens.dk.se
data.china.w5.education$preffalse.ind<-data.china.w5.education$sens.dk.rate-data.china.w5.education$nonsens.dk.rate
data.china.w5.education$preffalse.ind.se<-sqrt((data.china.w5.education$nonsens.dk.se^2)+(data.china.w5.education$sens.dk.se^2))
data.china.w5.education$preffalse.ind.u95ci<-data.china.w5.education$preffalse.ind+1.96*data.china.w5.education$preffalse.ind.se
data.china.w5.education$preffalse.ind.l95ci<-data.china.w5.education$preffalse.ind-1.96*data.china.w5.education$preffalse.ind.se
data.china.w5.education$conf.gov.se<-data.china.w5.education$conf.gov.sd/sqrt(data.china.w5.education$conf.gov.ob)
data.china.w5.education$conf.gov.u95ci<-data.china.w5.education$conf.gov.mean+1.96*data.china.w5.education$conf.gov.se
data.china.w5.education$conf.gov.l95ci<-data.china.w5.education$conf.gov.mean-1.96*data.china.w5.education$conf.gov.se
ggplot(data.china.w5.education, aes(x=education, y=preffalse.ind)) + xlab("Education") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5)  +  geom_segment(aes(x = education, y = preffalse.ind.l95ci, xend = education, yend = preffalse.ind.u95ci), color="grey60", alpha=.5, data = data.china.w5.education) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + geom_hline(aes(yintercept=0.00), colour="grey40", linetype="dotted") 
data.china.w5.education$variable<-"Education"
colnames(data.china.w5.education)[1] <- "level"

data.china.w5.institution<-aggregate(wvs.china.w5$conf.gov.ob, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)
data.china.w5.institution[,3]<-aggregate(wvs.china.w5$conf.gov.dk, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,4]<-aggregate(wvs.china.w5$resp.hr.ob, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,5]<-aggregate(wvs.china.w5$resp.hr.dk, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,6]<-aggregate(wvs.china.w5$democraticness.ob, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,7]<-aggregate(wvs.china.w5$democraticness.dk, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,8]<-aggregate(wvs.china.w5$conf.tv.ob, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,9]<-aggregate(wvs.china.w5$conf.tv.dk, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,10]<-aggregate(wvs.china.w5$gentrust.ob, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,11]<-aggregate(wvs.china.w5$gentrust.dk, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,12]<-aggregate(wvs.china.w5$lifesat.ob, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,13]<-aggregate(wvs.china.w5$lifesat.dk, by=list(wvs.china.w5$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w5.institution[,14]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$institution), FUN="mean", na.rm=TRUE)[2]
data.china.w5.institution[,15]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$institution), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w5.institution)<-c("institution","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w5.institution$conf.gov.dk.rate<-data.china.w5.institution$conf.gov.dk/data.china.w5.institution$conf.gov.ob
data.china.w5.institution$resp.hr.dk.rate<-data.china.w5.institution$resp.hr.dk/data.china.w5.institution$resp.hr.ob
data.china.w5.institution$democraticness.dk.rate<-data.china.w5.institution$democraticness.dk/data.china.w5.institution$democraticness.ob
data.china.w5.institution$conf.tv.dk.rate<-data.china.w5.institution$conf.tv.dk/data.china.w5.institution$conf.tv.ob
data.china.w5.institution$gentrust.dk.rate<-data.china.w5.institution$gentrust.dk/data.china.w5.institution$gentrust.ob
data.china.w5.institution$lifesat.dk.rate<-data.china.w5.institution$lifesat.dk/data.china.w5.institution$lifesat.ob
data.china.w5.institution$sens.ob<-data.china.w5.institution$conf.gov.ob+data.china.w5.institution$resp.hr.ob+data.china.w5.institution$democraticness.ob
data.china.w5.institution$sens.ob[data.china.w5.institution$conf.gov.ob==0]<-NA
data.china.w5.institution$sens.ob[data.china.w5.institution$resp.hr.ob==0]<-NA
data.china.w5.institution$sens.ob[data.china.w5.institution$democraticness.ob==0]<-NA
data.china.w5.institution$sens.dk<-data.china.w5.institution$conf.gov.dk+data.china.w5.institution$resp.hr.dk+data.china.w5.institution$democraticness.dk
data.china.w5.institution$sens.dk[data.china.w5.institution$conf.gov.ob==0]<-NA
data.china.w5.institution$sens.dk[data.china.w5.institution$resp.hr.ob==0]<-NA
data.china.w5.institution$sens.dk[data.china.w5.institution$democraticness.ob==0]<-NA
data.china.w5.institution$sens.dk.rate<-data.china.w5.institution$sens.dk/data.china.w5.institution$sens.ob
data.china.w5.institution$sens.dk.se<-sqrt((data.china.w5.institution$sens.dk.rate*(1-data.china.w5.institution$sens.dk.rate))/data.china.w5.institution$sens.ob)
data.china.w5.institution$sens.dk.u95ci<-data.china.w5.institution$sens.dk.rate+1.96*data.china.w5.institution$sens.dk.se
data.china.w5.institution$sens.dk.l95ci<-data.china.w5.institution$sens.dk.rate-1.96*data.china.w5.institution$sens.dk.se
data.china.w5.institution$nonsens.ob<-data.china.w5.institution$conf.tv.ob+data.china.w5.institution$lifesat.ob+data.china.w5.institution$gentrust.ob
data.china.w5.institution$nonsens.ob[data.china.w5.institution$conf.tv.ob==0]<-NA
data.china.w5.institution$nonsens.ob[data.china.w5.institution$lifesat.ob==0]<-NA
data.china.w5.institution$nonsens.ob[data.china.w5.institution$gentrust.ob==0]<-NA
data.china.w5.institution$nonsens.dk<-data.china.w5.institution$conf.tv.dk+data.china.w5.institution$lifesat.dk+data.china.w5.institution$gentrust.dk
data.china.w5.institution$nonsens.dk[data.china.w5.institution$conf.tv.ob==0]<-NA
data.china.w5.institution$nonsens.dk[data.china.w5.institution$lifesat.ob==0]<-NA
data.china.w5.institution$nonsens.dk[data.china.w5.institution$gentrust.ob==0]<-NA
data.china.w5.institution$nonsens.dk.rate<-data.china.w5.institution$nonsens.dk/data.china.w5.institution$nonsens.ob
data.china.w5.institution$nonsens.dk.se<-sqrt((data.china.w5.institution$nonsens.dk.rate*(1-data.china.w5.institution$nonsens.dk.rate))/data.china.w5.institution$nonsens.ob)
data.china.w5.institution$nonsens.dk.u95ci<-data.china.w5.institution$nonsens.dk.rate+1.96*data.china.w5.institution$nonsens.dk.se
data.china.w5.institution$nonsens.dk.l95ci<-data.china.w5.institution$nonsens.dk.rate-1.96*data.china.w5.institution$nonsens.dk.se
data.china.w5.institution$preffalse.ind<-data.china.w5.institution$sens.dk.rate-data.china.w5.institution$nonsens.dk.rate
data.china.w5.institution$preffalse.ind.se<-sqrt((data.china.w5.institution$nonsens.dk.se^2)+(data.china.w5.institution$sens.dk.se^2))
data.china.w5.institution$preffalse.ind.u95ci<-data.china.w5.institution$preffalse.ind+1.96*data.china.w5.institution$preffalse.ind.se
data.china.w5.institution$preffalse.ind.l95ci<-data.china.w5.institution$preffalse.ind-1.96*data.china.w5.institution$preffalse.ind.se
data.china.w5.institution$conf.gov.se<-data.china.w5.institution$conf.gov.sd/sqrt(data.china.w5.institution$conf.gov.ob)
data.china.w5.institution$conf.gov.u95ci<-data.china.w5.institution$conf.gov.mean+1.96*data.china.w5.institution$conf.gov.se
data.china.w5.institution$conf.gov.l95ci<-data.china.w5.institution$conf.gov.mean-1.96*data.china.w5.institution$conf.gov.se
ggplot(data.china.w5.institution, aes(x=institution, y=preffalse.ind)) + xlab("institution") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5)  +  geom_segment(aes(x = institution, y = preffalse.ind.l95ci, xend = institution, yend = preffalse.ind.u95ci), color="grey60", alpha=.5, data = data.china.w5.institution) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + geom_hline(aes(yintercept=0.00), colour="grey40", linetype="dotted") 
data.china.w5.institution$variable<-"Institution"
colnames(data.china.w5.institution)[1] <- "level"

data.china.w5.class<-aggregate(wvs.china.w5$conf.gov.ob, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)
data.china.w5.class[,3]<-aggregate(wvs.china.w5$conf.gov.dk, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,4]<-aggregate(wvs.china.w5$resp.hr.ob, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,5]<-aggregate(wvs.china.w5$resp.hr.dk, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,6]<-aggregate(wvs.china.w5$democraticness.ob, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,7]<-aggregate(wvs.china.w5$democraticness.dk, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,8]<-aggregate(wvs.china.w5$conf.tv.ob, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,9]<-aggregate(wvs.china.w5$conf.tv.dk, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,10]<-aggregate(wvs.china.w5$gentrust.ob, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,11]<-aggregate(wvs.china.w5$gentrust.dk, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,12]<-aggregate(wvs.china.w5$lifesat.ob, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,13]<-aggregate(wvs.china.w5$lifesat.dk, by=list(wvs.china.w5$class), FUN="sum", na.rm=TRUE)[2]
data.china.w5.class[,14]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$class), FUN="mean", na.rm=TRUE)[2]
data.china.w5.class[,15]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$class), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w5.class)<-c("class","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w5.class$conf.gov.dk.rate<-data.china.w5.class$conf.gov.dk/data.china.w5.class$conf.gov.ob
data.china.w5.class$resp.hr.dk.rate<-data.china.w5.class$resp.hr.dk/data.china.w5.class$resp.hr.ob
data.china.w5.class$democraticness.dk.rate<-data.china.w5.class$democraticness.dk/data.china.w5.class$democraticness.ob
data.china.w5.class$conf.tv.dk.rate<-data.china.w5.class$conf.tv.dk/data.china.w5.class$conf.tv.ob
data.china.w5.class$gentrust.dk.rate<-data.china.w5.class$gentrust.dk/data.china.w5.class$gentrust.ob
data.china.w5.class$lifesat.dk.rate<-data.china.w5.class$lifesat.dk/data.china.w5.class$lifesat.ob
data.china.w5.class$sens.ob<-data.china.w5.class$conf.gov.ob+data.china.w5.class$resp.hr.ob+data.china.w5.class$democraticness.ob
data.china.w5.class$sens.ob[data.china.w5.class$conf.gov.ob==0]<-NA
data.china.w5.class$sens.ob[data.china.w5.class$resp.hr.ob==0]<-NA
data.china.w5.class$sens.ob[data.china.w5.class$democraticness.ob==0]<-NA
data.china.w5.class$sens.dk<-data.china.w5.class$conf.gov.dk+data.china.w5.class$resp.hr.dk+data.china.w5.class$democraticness.dk
data.china.w5.class$sens.dk[data.china.w5.class$conf.gov.ob==0]<-NA
data.china.w5.class$sens.dk[data.china.w5.class$resp.hr.ob==0]<-NA
data.china.w5.class$sens.dk[data.china.w5.class$democraticness.ob==0]<-NA
data.china.w5.class$sens.dk.rate<-data.china.w5.class$sens.dk/data.china.w5.class$sens.ob
data.china.w5.class$sens.dk.se<-sqrt((data.china.w5.class$sens.dk.rate*(1-data.china.w5.class$sens.dk.rate))/data.china.w5.class$sens.ob)
data.china.w5.class$sens.dk.u95ci<-data.china.w5.class$sens.dk.rate+1.96*data.china.w5.class$sens.dk.se
data.china.w5.class$sens.dk.l95ci<-data.china.w5.class$sens.dk.rate-1.96*data.china.w5.class$sens.dk.se
data.china.w5.class$nonsens.ob<-data.china.w5.class$conf.tv.ob+data.china.w5.class$lifesat.ob+data.china.w5.class$gentrust.ob
data.china.w5.class$nonsens.ob[data.china.w5.class$conf.tv.ob==0]<-NA
data.china.w5.class$nonsens.ob[data.china.w5.class$lifesat.ob==0]<-NA
data.china.w5.class$nonsens.ob[data.china.w5.class$gentrust.ob==0]<-NA
data.china.w5.class$nonsens.dk<-data.china.w5.class$conf.tv.dk+data.china.w5.class$lifesat.dk+data.china.w5.class$gentrust.dk
data.china.w5.class$nonsens.dk[data.china.w5.class$conf.tv.ob==0]<-NA
data.china.w5.class$nonsens.dk[data.china.w5.class$lifesat.ob==0]<-NA
data.china.w5.class$nonsens.dk[data.china.w5.class$gentrust.ob==0]<-NA
data.china.w5.class$nonsens.dk.rate<-data.china.w5.class$nonsens.dk/data.china.w5.class$nonsens.ob
data.china.w5.class$nonsens.dk.se<-sqrt((data.china.w5.class$nonsens.dk.rate*(1-data.china.w5.class$nonsens.dk.rate))/data.china.w5.class$nonsens.ob)
data.china.w5.class$nonsens.dk.u95ci<-data.china.w5.class$nonsens.dk.rate+1.96*data.china.w5.class$nonsens.dk.se
data.china.w5.class$nonsens.dk.l95ci<-data.china.w5.class$nonsens.dk.rate-1.96*data.china.w5.class$nonsens.dk.se
data.china.w5.class$preffalse.ind<-data.china.w5.class$sens.dk.rate-data.china.w5.class$nonsens.dk.rate
data.china.w5.class$preffalse.ind.se<-sqrt((data.china.w5.class$nonsens.dk.se^2)+(data.china.w5.class$sens.dk.se^2))
data.china.w5.class$preffalse.ind.u95ci<-data.china.w5.class$preffalse.ind+1.96*data.china.w5.class$preffalse.ind.se
data.china.w5.class$preffalse.ind.l95ci<-data.china.w5.class$preffalse.ind-1.96*data.china.w5.class$preffalse.ind.se
data.china.w5.class$conf.gov.se<-data.china.w5.class$conf.gov.sd/sqrt(data.china.w5.class$conf.gov.ob)
data.china.w5.class$conf.gov.u95ci<-data.china.w5.class$conf.gov.mean+1.96*data.china.w5.class$conf.gov.se
data.china.w5.class$conf.gov.l95ci<-data.china.w5.class$conf.gov.mean-1.96*data.china.w5.class$conf.gov.se
ggplot(data.china.w5.class, aes(x=class, y=preffalse.ind)) + xlab("class") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5)  +  geom_segment(aes(x = class, y = preffalse.ind.l95ci, xend = class, yend = preffalse.ind.u95ci), color="grey60", alpha=.5, data = data.china.w5.class) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + geom_hline(aes(yintercept=0.00), colour="grey40", linetype="dotted") 
data.china.w5.class$variable<-"Class"
colnames(data.china.w5.class)[1] <- "level"

data.china.w5<-rbind(data.china.w5.education, data.china.w5.gender, data.china.w5.age, data.china.w5.class, data.china.w5.institution)


#WVS - Wave 6
wvs.china<-subset(wvs, wvs$S003=="China")
wvs.china.w6<-subset(wvs.china, wvs.china$S025=="China (2012)")

wvs.china.w6$E069_11
wvs.china.w6$conf.gov<-NA
wvs.china.w6$conf.gov[wvs.china.w6$E069_11=="A great deal"]<-4
wvs.china.w6$conf.gov[wvs.china.w6$E069_11=="Quite a lot"]<-3
wvs.china.w6$conf.gov[wvs.china.w6$E069_11=="Not very much"]<-2
wvs.china.w6$conf.gov[wvs.china.w6$E069_11=="None at all"]<-1

wvs.china.w6$birthyear<-NA
wvs.china.w6$birthyear[as.numeric(paste(wvs.china.w6$X002))<2000]<-"1980+"
wvs.china.w6$birthyear[as.numeric(paste(wvs.china.w6$X002))<1980]<-"1970s"
wvs.china.w6$birthyear[as.numeric(paste(wvs.china.w6$X002))<1970]<-"1960s"
wvs.china.w6$birthyear[as.numeric(paste(wvs.china.w6$X002))<1960]<-"1950s"
wvs.china.w6$birthyear[as.numeric(paste(wvs.china.w6$X002))<1950]<-"<1950"

wvs.china.w6$institution<-NA
wvs.china.w6$institution[wvs.china.w6$X052=="Public institution"]<-"Public"
wvs.china.w6$institution[wvs.china.w6$X052=="Self-employed"]<-"Private"
wvs.china.w6$institution[wvs.china.w6$X052=="Private business"]<-"Private"
wvs.china.w6$institution[wvs.china.w6$X052=="Private non-profit organization"]<-"Private"

wvs.china.w6$education<-NA
wvs.china.w6$education[wvs.china.w6$X025=="Complete secondary school: technical/vocational type/Seconda"]<-"Secondary"
wvs.china.w6$education[wvs.china.w6$X025=="Complete secondary: university-preparatory type/Full seconda"]<-"Secondary"
wvs.china.w6$education[wvs.china.w6$X025=="Completed (compulsory) elementary education"]<-"Elementary"
wvs.china.w6$education[wvs.china.w6$X025=="University with degree/Higher education - upper-level tertia"]<-"University"
wvs.china.w6$education[wvs.china.w6$X025=="Not applicable; No formal education"]<-"No education"

wvs.china.w6$class<-NA
wvs.china.w6$class[wvs.china.w6$X045=="Lower class"]<-"1 - Lower class"
wvs.china.w6$class[wvs.china.w6$X045=="Working class"]<-"2 - Lower middle class"
wvs.china.w6$class[wvs.china.w6$X045=="Lower middle class"]<-"3 - Middle class"
wvs.china.w6$class[wvs.china.w6$X045=="Upper middle class"]<-"4 - Upper middle class"
wvs.china.w6$class[wvs.china.w6$X045=="Upper class"]<-"5 - Upper class"

data.china.w6.age<-aggregate(wvs.china.w6$conf.gov.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)
data.china.w6.age[,3]<-aggregate(wvs.china.w6$conf.gov.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,4]<-aggregate(wvs.china.w6$resp.hr.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,5]<-aggregate(wvs.china.w6$resp.hr.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,6]<-aggregate(wvs.china.w6$democraticness.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,7]<-aggregate(wvs.china.w6$democraticness.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,8]<-aggregate(wvs.china.w6$conf.tv.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,9]<-aggregate(wvs.china.w6$conf.tv.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,10]<-aggregate(wvs.china.w6$gentrust.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,11]<-aggregate(wvs.china.w6$gentrust.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,12]<-aggregate(wvs.china.w6$lifesat.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,13]<-aggregate(wvs.china.w6$lifesat.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,14]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$birthyear), FUN="mean", na.rm=TRUE)[2]
data.china.w6.age[,15]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$birthyear), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w6.age)<-c("birthyear","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w6.age$conf.gov.dk.rate<-data.china.w6.age$conf.gov.dk/data.china.w6.age$conf.gov.ob
data.china.w6.age$resp.hr.dk.rate<-data.china.w6.age$resp.hr.dk/data.china.w6.age$resp.hr.ob
data.china.w6.age$democraticness.dk.rate<-data.china.w6.age$democraticness.dk/data.china.w6.age$democraticness.ob
data.china.w6.age$conf.tv.dk.rate<-data.china.w6.age$conf.tv.dk/data.china.w6.age$conf.tv.ob
data.china.w6.age$gentrust.dk.rate<-data.china.w6.age$gentrust.dk/data.china.w6.age$gentrust.ob
data.china.w6.age$lifesat.dk.rate<-data.china.w6.age$lifesat.dk/data.china.w6.age$lifesat.ob
data.china.w6.age$sens.ob<-data.china.w6.age$conf.gov.ob+data.china.w6.age$resp.hr.ob+data.china.w6.age$democraticness.ob
data.china.w6.age$sens.ob[data.china.w6.age$conf.gov.ob==0]<-NA
data.china.w6.age$sens.ob[data.china.w6.age$resp.hr.ob==0]<-NA
data.china.w6.age$sens.ob[data.china.w6.age$democraticness.ob==0]<-NA
data.china.w6.age$sens.dk<-data.china.w6.age$conf.gov.dk+data.china.w6.age$resp.hr.dk+data.china.w6.age$democraticness.dk
data.china.w6.age$sens.dk[data.china.w6.age$conf.gov.ob==0]<-NA
data.china.w6.age$sens.dk[data.china.w6.age$resp.hr.ob==0]<-NA
data.china.w6.age$sens.dk[data.china.w6.age$democraticness.ob==0]<-NA
data.china.w6.age$sens.dk.rate<-data.china.w6.age$sens.dk/data.china.w6.age$sens.ob
data.china.w6.age$sens.dk.se<-sqrt((data.china.w6.age$sens.dk.rate*(1-data.china.w6.age$sens.dk.rate))/data.china.w6.age$sens.ob)
data.china.w6.age$sens.dk.u95ci<-data.china.w6.age$sens.dk.rate+1.96*data.china.w6.age$sens.dk.se
data.china.w6.age$sens.dk.l95ci<-data.china.w6.age$sens.dk.rate-1.96*data.china.w6.age$sens.dk.se
data.china.w6.age$nonsens.ob<-data.china.w6.age$conf.tv.ob+data.china.w6.age$lifesat.ob+data.china.w6.age$gentrust.ob
data.china.w6.age$nonsens.ob[data.china.w6.age$conf.tv.ob==0]<-NA
data.china.w6.age$nonsens.ob[data.china.w6.age$lifesat.ob==0]<-NA
data.china.w6.age$nonsens.ob[data.china.w6.age$gentrust.ob==0]<-NA
data.china.w6.age$nonsens.dk<-data.china.w6.age$conf.tv.dk+data.china.w6.age$lifesat.dk+data.china.w6.age$gentrust.dk
data.china.w6.age$nonsens.dk[data.china.w6.age$conf.tv.ob==0]<-NA
data.china.w6.age$nonsens.dk[data.china.w6.age$lifesat.ob==0]<-NA
data.china.w6.age$nonsens.dk[data.china.w6.age$gentrust.ob==0]<-NA
data.china.w6.age$nonsens.dk.rate<-data.china.w6.age$nonsens.dk/data.china.w6.age$nonsens.ob
data.china.w6.age$nonsens.dk.se<-sqrt((data.china.w6.age$nonsens.dk.rate*(1-data.china.w6.age$nonsens.dk.rate))/data.china.w6.age$nonsens.ob)
data.china.w6.age$nonsens.dk.u95ci<-data.china.w6.age$nonsens.dk.rate+1.96*data.china.w6.age$nonsens.dk.se
data.china.w6.age$nonsens.dk.l95ci<-data.china.w6.age$nonsens.dk.rate-1.96*data.china.w6.age$nonsens.dk.se
data.china.w6.age$preffalse.ind<-data.china.w6.age$sens.dk.rate-data.china.w6.age$nonsens.dk.rate
data.china.w6.age$preffalse.ind.se<-sqrt((data.china.w6.age$nonsens.dk.se^2)+(data.china.w6.age$sens.dk.se^2))
data.china.w6.age$preffalse.ind.u95ci<-data.china.w6.age$preffalse.ind+1.96*data.china.w6.age$preffalse.ind.se
data.china.w6.age$preffalse.ind.l95ci<-data.china.w6.age$preffalse.ind-1.96*data.china.w6.age$preffalse.ind.se
data.china.w6.age$conf.gov.se<-data.china.w6.age$conf.gov.sd/sqrt(data.china.w6.age$conf.gov.ob)
data.china.w6.age$conf.gov.u95ci<-data.china.w6.age$conf.gov.mean+1.96*data.china.w6.age$conf.gov.se
data.china.w6.age$conf.gov.l95ci<-data.china.w6.age$conf.gov.mean-1.96*data.china.w6.age$conf.gov.se
data.china.w6.age$variable<-"Birthyear"
colnames(data.china.w6.age)[1] <- "level"

data.china.w6.gender<-aggregate(wvs.china.w6$conf.gov.ob, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)
data.china.w6.gender[,3]<-aggregate(wvs.china.w6$conf.gov.dk, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,4]<-aggregate(wvs.china.w6$resp.hr.ob, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,5]<-aggregate(wvs.china.w6$resp.hr.dk, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,6]<-aggregate(wvs.china.w6$democraticness.ob, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,7]<-aggregate(wvs.china.w6$democraticness.dk, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,8]<-aggregate(wvs.china.w6$conf.tv.ob, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,9]<-aggregate(wvs.china.w6$conf.tv.dk, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,10]<-aggregate(wvs.china.w6$gentrust.ob, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,11]<-aggregate(wvs.china.w6$gentrust.dk, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,12]<-aggregate(wvs.china.w6$lifesat.ob, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,13]<-aggregate(wvs.china.w6$lifesat.dk, by=list(wvs.china.w6$X001), FUN="sum", na.rm=TRUE)[2]
data.china.w6.gender[,14]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$X001), FUN="mean", na.rm=TRUE)[2]
data.china.w6.gender[,15]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$X001), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w6.gender)<-c("gender","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w6.gender$conf.gov.dk.rate<-data.china.w6.gender$conf.gov.dk/data.china.w6.gender$conf.gov.ob
data.china.w6.gender$resp.hr.dk.rate<-data.china.w6.gender$resp.hr.dk/data.china.w6.gender$resp.hr.ob
data.china.w6.gender$democraticness.dk.rate<-data.china.w6.gender$democraticness.dk/data.china.w6.gender$democraticness.ob
data.china.w6.gender$conf.tv.dk.rate<-data.china.w6.gender$conf.tv.dk/data.china.w6.gender$conf.tv.ob
data.china.w6.gender$gentrust.dk.rate<-data.china.w6.gender$gentrust.dk/data.china.w6.gender$gentrust.ob
data.china.w6.gender$lifesat.dk.rate<-data.china.w6.gender$lifesat.dk/data.china.w6.gender$lifesat.ob
data.china.w6.gender$sens.ob<-data.china.w6.gender$conf.gov.ob+data.china.w6.gender$resp.hr.ob+data.china.w6.gender$democraticness.ob
data.china.w6.gender$sens.ob[data.china.w6.gender$conf.gov.ob==0]<-NA
data.china.w6.gender$sens.ob[data.china.w6.gender$resp.hr.ob==0]<-NA
data.china.w6.gender$sens.ob[data.china.w6.gender$democraticness.ob==0]<-NA
data.china.w6.gender$sens.dk<-data.china.w6.gender$conf.gov.dk+data.china.w6.gender$resp.hr.dk+data.china.w6.gender$democraticness.dk
data.china.w6.gender$sens.dk[data.china.w6.gender$conf.gov.ob==0]<-NA
data.china.w6.gender$sens.dk[data.china.w6.gender$resp.hr.ob==0]<-NA
data.china.w6.gender$sens.dk[data.china.w6.gender$democraticness.ob==0]<-NA
data.china.w6.gender$sens.dk.rate<-data.china.w6.gender$sens.dk/data.china.w6.gender$sens.ob
data.china.w6.gender$sens.dk.se<-sqrt((data.china.w6.gender$sens.dk.rate*(1-data.china.w6.gender$sens.dk.rate))/data.china.w6.gender$sens.ob)
data.china.w6.gender$sens.dk.u95ci<-data.china.w6.gender$sens.dk.rate+1.96*data.china.w6.gender$sens.dk.se
data.china.w6.gender$sens.dk.l95ci<-data.china.w6.gender$sens.dk.rate-1.96*data.china.w6.gender$sens.dk.se
data.china.w6.gender$nonsens.ob<-data.china.w6.gender$conf.tv.ob+data.china.w6.gender$lifesat.ob+data.china.w6.gender$gentrust.ob
data.china.w6.gender$nonsens.ob[data.china.w6.gender$conf.tv.ob==0]<-NA
data.china.w6.gender$nonsens.ob[data.china.w6.gender$lifesat.ob==0]<-NA
data.china.w6.gender$nonsens.ob[data.china.w6.gender$gentrust.ob==0]<-NA
data.china.w6.gender$nonsens.dk<-data.china.w6.gender$conf.tv.dk+data.china.w6.gender$lifesat.dk+data.china.w6.gender$gentrust.dk
data.china.w6.gender$nonsens.dk[data.china.w6.gender$conf.tv.ob==0]<-NA
data.china.w6.gender$nonsens.dk[data.china.w6.gender$lifesat.ob==0]<-NA
data.china.w6.gender$nonsens.dk[data.china.w6.gender$gentrust.ob==0]<-NA
data.china.w6.gender$nonsens.dk.rate<-data.china.w6.gender$nonsens.dk/data.china.w6.gender$nonsens.ob
data.china.w6.gender$nonsens.dk.se<-sqrt((data.china.w6.gender$nonsens.dk.rate*(1-data.china.w6.gender$nonsens.dk.rate))/data.china.w6.gender$nonsens.ob)
data.china.w6.gender$nonsens.dk.u95ci<-data.china.w6.gender$nonsens.dk.rate+1.96*data.china.w6.gender$nonsens.dk.se
data.china.w6.gender$nonsens.dk.l95ci<-data.china.w6.gender$nonsens.dk.rate-1.96*data.china.w6.gender$nonsens.dk.se
data.china.w6.gender$preffalse.ind<-data.china.w6.gender$sens.dk.rate-data.china.w6.gender$nonsens.dk.rate
data.china.w6.gender$preffalse.ind.se<-sqrt((data.china.w6.gender$nonsens.dk.se^2)+(data.china.w6.gender$sens.dk.se^2))
data.china.w6.gender$preffalse.ind.u95ci<-data.china.w6.gender$preffalse.ind+1.96*data.china.w6.gender$preffalse.ind.se
data.china.w6.gender$preffalse.ind.l95ci<-data.china.w6.gender$preffalse.ind-1.96*data.china.w6.gender$preffalse.ind.se
data.china.w6.gender$conf.gov.se<-data.china.w6.gender$conf.gov.sd/sqrt(data.china.w6.gender$conf.gov.ob)
data.china.w6.gender$conf.gov.u95ci<-data.china.w6.gender$conf.gov.mean+1.96*data.china.w6.gender$conf.gov.se
data.china.w6.gender$conf.gov.l95ci<-data.china.w6.gender$conf.gov.mean-1.96*data.china.w6.gender$conf.gov.se
ggplot(data.china.w6.gender, aes(x=gender, y=preffalse.ind)) + xlab("Gender") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5)  +  geom_segment(aes(x = gender, y = preffalse.ind.l95ci, xend = gender, yend = preffalse.ind.u95ci), color="grey60", alpha=.5, data = data.china.w6.gender) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + geom_hline(aes(yintercept=0.00), colour="grey40", linetype="dotted") 
data.china.w6.gender$variable<-"Gender"
colnames(data.china.w6.gender)[1] <- "level"

data.china.w6.education<-aggregate(wvs.china.w6$conf.gov.ob, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)
data.china.w6.education[,3]<-aggregate(wvs.china.w6$conf.gov.dk, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,4]<-aggregate(wvs.china.w6$resp.hr.ob, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,5]<-aggregate(wvs.china.w6$resp.hr.dk, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,6]<-aggregate(wvs.china.w6$democraticness.ob, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,7]<-aggregate(wvs.china.w6$democraticness.dk, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,8]<-aggregate(wvs.china.w6$conf.tv.ob, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,9]<-aggregate(wvs.china.w6$conf.tv.dk, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,10]<-aggregate(wvs.china.w6$gentrust.ob, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,11]<-aggregate(wvs.china.w6$gentrust.dk, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,12]<-aggregate(wvs.china.w6$lifesat.ob, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,13]<-aggregate(wvs.china.w6$lifesat.dk, by=list(wvs.china.w6$education), FUN="sum", na.rm=TRUE)[2]
data.china.w6.education[,14]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$education), FUN="mean", na.rm=TRUE)[2]
data.china.w6.education[,15]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$education), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w6.education)<-c("education","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w6.education$conf.gov.dk.rate<-data.china.w6.education$conf.gov.dk/data.china.w6.education$conf.gov.ob
data.china.w6.education$resp.hr.dk.rate<-data.china.w6.education$resp.hr.dk/data.china.w6.education$resp.hr.ob
data.china.w6.education$democraticness.dk.rate<-data.china.w6.education$democraticness.dk/data.china.w6.education$democraticness.ob
data.china.w6.education$conf.tv.dk.rate<-data.china.w6.education$conf.tv.dk/data.china.w6.education$conf.tv.ob
data.china.w6.education$gentrust.dk.rate<-data.china.w6.education$gentrust.dk/data.china.w6.education$gentrust.ob
data.china.w6.education$lifesat.dk.rate<-data.china.w6.education$lifesat.dk/data.china.w6.education$lifesat.ob
data.china.w6.education$sens.ob<-data.china.w6.education$conf.gov.ob+data.china.w6.education$resp.hr.ob+data.china.w6.education$democraticness.ob
data.china.w6.education$sens.ob[data.china.w6.education$conf.gov.ob==0]<-NA
data.china.w6.education$sens.ob[data.china.w6.education$resp.hr.ob==0]<-NA
data.china.w6.education$sens.ob[data.china.w6.education$democraticness.ob==0]<-NA
data.china.w6.education$sens.dk<-data.china.w6.education$conf.gov.dk+data.china.w6.education$resp.hr.dk+data.china.w6.education$democraticness.dk
data.china.w6.education$sens.dk[data.china.w6.education$conf.gov.ob==0]<-NA
data.china.w6.education$sens.dk[data.china.w6.education$resp.hr.ob==0]<-NA
data.china.w6.education$sens.dk[data.china.w6.education$democraticness.ob==0]<-NA
data.china.w6.education$sens.dk.rate<-data.china.w6.education$sens.dk/data.china.w6.education$sens.ob
data.china.w6.education$sens.dk.se<-sqrt((data.china.w6.education$sens.dk.rate*(1-data.china.w6.education$sens.dk.rate))/data.china.w6.education$sens.ob)
data.china.w6.education$sens.dk.u95ci<-data.china.w6.education$sens.dk.rate+1.96*data.china.w6.education$sens.dk.se
data.china.w6.education$sens.dk.l95ci<-data.china.w6.education$sens.dk.rate-1.96*data.china.w6.education$sens.dk.se
data.china.w6.education$nonsens.ob<-data.china.w6.education$conf.tv.ob+data.china.w6.education$lifesat.ob+data.china.w6.education$gentrust.ob
data.china.w6.education$nonsens.ob[data.china.w6.education$conf.tv.ob==0]<-NA
data.china.w6.education$nonsens.ob[data.china.w6.education$lifesat.ob==0]<-NA
data.china.w6.education$nonsens.ob[data.china.w6.education$gentrust.ob==0]<-NA
data.china.w6.education$nonsens.dk<-data.china.w6.education$conf.tv.dk+data.china.w6.education$lifesat.dk+data.china.w6.education$gentrust.dk
data.china.w6.education$nonsens.dk[data.china.w6.education$conf.tv.ob==0]<-NA
data.china.w6.education$nonsens.dk[data.china.w6.education$lifesat.ob==0]<-NA
data.china.w6.education$nonsens.dk[data.china.w6.education$gentrust.ob==0]<-NA
data.china.w6.education$nonsens.dk.rate<-data.china.w6.education$nonsens.dk/data.china.w6.education$nonsens.ob
data.china.w6.education$nonsens.dk.se<-sqrt((data.china.w6.education$nonsens.dk.rate*(1-data.china.w6.education$nonsens.dk.rate))/data.china.w6.education$nonsens.ob)
data.china.w6.education$nonsens.dk.u95ci<-data.china.w6.education$nonsens.dk.rate+1.96*data.china.w6.education$nonsens.dk.se
data.china.w6.education$nonsens.dk.l95ci<-data.china.w6.education$nonsens.dk.rate-1.96*data.china.w6.education$nonsens.dk.se
data.china.w6.education$preffalse.ind<-data.china.w6.education$sens.dk.rate-data.china.w6.education$nonsens.dk.rate
data.china.w6.education$preffalse.ind.se<-sqrt((data.china.w6.education$nonsens.dk.se^2)+(data.china.w6.education$sens.dk.se^2))
data.china.w6.education$preffalse.ind.u95ci<-data.china.w6.education$preffalse.ind+1.96*data.china.w6.education$preffalse.ind.se
data.china.w6.education$preffalse.ind.l95ci<-data.china.w6.education$preffalse.ind-1.96*data.china.w6.education$preffalse.ind.se
data.china.w6.education$conf.gov.se<-data.china.w6.education$conf.gov.sd/sqrt(data.china.w6.education$conf.gov.ob)
data.china.w6.education$conf.gov.u95ci<-data.china.w6.education$conf.gov.mean+1.96*data.china.w6.education$conf.gov.se
data.china.w6.education$conf.gov.l95ci<-data.china.w6.education$conf.gov.mean-1.96*data.china.w6.education$conf.gov.se
ggplot(data.china.w6.education, aes(x=education, y=preffalse.ind)) + xlab("Education") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5)  +  geom_segment(aes(x = education, y = preffalse.ind.l95ci, xend = education, yend = preffalse.ind.u95ci), color="grey60", alpha=.5, data = data.china.w6.education) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + geom_hline(aes(yintercept=0.00), colour="grey40", linetype="dotted") 
data.china.w6.education$variable<-"Education"
colnames(data.china.w6.education)[1] <- "level"

data.china.w6.institution<-aggregate(wvs.china.w6$conf.gov.ob, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)
data.china.w6.institution[,3]<-aggregate(wvs.china.w6$conf.gov.dk, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,4]<-aggregate(wvs.china.w6$resp.hr.ob, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,5]<-aggregate(wvs.china.w6$resp.hr.dk, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,6]<-aggregate(wvs.china.w6$democraticness.ob, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,7]<-aggregate(wvs.china.w6$democraticness.dk, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,8]<-aggregate(wvs.china.w6$conf.tv.ob, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,9]<-aggregate(wvs.china.w6$conf.tv.dk, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,10]<-aggregate(wvs.china.w6$gentrust.ob, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,11]<-aggregate(wvs.china.w6$gentrust.dk, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,12]<-aggregate(wvs.china.w6$lifesat.ob, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,13]<-aggregate(wvs.china.w6$lifesat.dk, by=list(wvs.china.w6$institution), FUN="sum", na.rm=TRUE)[2]
data.china.w6.institution[,14]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$institution), FUN="mean", na.rm=TRUE)[2]
data.china.w6.institution[,15]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$institution), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w6.institution)<-c("institution","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w6.institution$conf.gov.dk.rate<-data.china.w6.institution$conf.gov.dk/data.china.w6.institution$conf.gov.ob
data.china.w6.institution$resp.hr.dk.rate<-data.china.w6.institution$resp.hr.dk/data.china.w6.institution$resp.hr.ob
data.china.w6.institution$democraticness.dk.rate<-data.china.w6.institution$democraticness.dk/data.china.w6.institution$democraticness.ob
data.china.w6.institution$conf.tv.dk.rate<-data.china.w6.institution$conf.tv.dk/data.china.w6.institution$conf.tv.ob
data.china.w6.institution$gentrust.dk.rate<-data.china.w6.institution$gentrust.dk/data.china.w6.institution$gentrust.ob
data.china.w6.institution$lifesat.dk.rate<-data.china.w6.institution$lifesat.dk/data.china.w6.institution$lifesat.ob
data.china.w6.institution$sens.ob<-data.china.w6.institution$conf.gov.ob+data.china.w6.institution$resp.hr.ob+data.china.w6.institution$democraticness.ob
data.china.w6.institution$sens.ob[data.china.w6.institution$conf.gov.ob==0]<-NA
data.china.w6.institution$sens.ob[data.china.w6.institution$resp.hr.ob==0]<-NA
data.china.w6.institution$sens.ob[data.china.w6.institution$democraticness.ob==0]<-NA
data.china.w6.institution$sens.dk<-data.china.w6.institution$conf.gov.dk+data.china.w6.institution$resp.hr.dk+data.china.w6.institution$democraticness.dk
data.china.w6.institution$sens.dk[data.china.w6.institution$conf.gov.ob==0]<-NA
data.china.w6.institution$sens.dk[data.china.w6.institution$resp.hr.ob==0]<-NA
data.china.w6.institution$sens.dk[data.china.w6.institution$democraticness.ob==0]<-NA
data.china.w6.institution$sens.dk.rate<-data.china.w6.institution$sens.dk/data.china.w6.institution$sens.ob
data.china.w6.institution$sens.dk.se<-sqrt((data.china.w6.institution$sens.dk.rate*(1-data.china.w6.institution$sens.dk.rate))/data.china.w6.institution$sens.ob)
data.china.w6.institution$sens.dk.u95ci<-data.china.w6.institution$sens.dk.rate+1.96*data.china.w6.institution$sens.dk.se
data.china.w6.institution$sens.dk.l95ci<-data.china.w6.institution$sens.dk.rate-1.96*data.china.w6.institution$sens.dk.se
data.china.w6.institution$nonsens.ob<-data.china.w6.institution$conf.tv.ob+data.china.w6.institution$lifesat.ob+data.china.w6.institution$gentrust.ob
data.china.w6.institution$nonsens.ob[data.china.w6.institution$conf.tv.ob==0]<-NA
data.china.w6.institution$nonsens.ob[data.china.w6.institution$lifesat.ob==0]<-NA
data.china.w6.institution$nonsens.ob[data.china.w6.institution$gentrust.ob==0]<-NA
data.china.w6.institution$nonsens.dk<-data.china.w6.institution$conf.tv.dk+data.china.w6.institution$lifesat.dk+data.china.w6.institution$gentrust.dk
data.china.w6.institution$nonsens.dk[data.china.w6.institution$conf.tv.ob==0]<-NA
data.china.w6.institution$nonsens.dk[data.china.w6.institution$lifesat.ob==0]<-NA
data.china.w6.institution$nonsens.dk[data.china.w6.institution$gentrust.ob==0]<-NA
data.china.w6.institution$nonsens.dk.rate<-data.china.w6.institution$nonsens.dk/data.china.w6.institution$nonsens.ob
data.china.w6.institution$nonsens.dk.se<-sqrt((data.china.w6.institution$nonsens.dk.rate*(1-data.china.w6.institution$nonsens.dk.rate))/data.china.w6.institution$nonsens.ob)
data.china.w6.institution$nonsens.dk.u95ci<-data.china.w6.institution$nonsens.dk.rate+1.96*data.china.w6.institution$nonsens.dk.se
data.china.w6.institution$nonsens.dk.l95ci<-data.china.w6.institution$nonsens.dk.rate-1.96*data.china.w6.institution$nonsens.dk.se
data.china.w6.institution$preffalse.ind<-data.china.w6.institution$sens.dk.rate-data.china.w6.institution$nonsens.dk.rate
data.china.w6.institution$preffalse.ind.se<-sqrt((data.china.w6.institution$nonsens.dk.se^2)+(data.china.w6.institution$sens.dk.se^2))
data.china.w6.institution$preffalse.ind.u95ci<-data.china.w6.institution$preffalse.ind+1.96*data.china.w6.institution$preffalse.ind.se
data.china.w6.institution$preffalse.ind.l95ci<-data.china.w6.institution$preffalse.ind-1.96*data.china.w6.institution$preffalse.ind.se
data.china.w6.institution$conf.gov.se<-data.china.w6.institution$conf.gov.sd/sqrt(data.china.w6.institution$conf.gov.ob)
data.china.w6.institution$conf.gov.u95ci<-data.china.w6.institution$conf.gov.mean+1.96*data.china.w6.institution$conf.gov.se
data.china.w6.institution$conf.gov.l95ci<-data.china.w6.institution$conf.gov.mean-1.96*data.china.w6.institution$conf.gov.se
ggplot(data.china.w6.institution, aes(x=institution, y=preffalse.ind)) + xlab("institution") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5)  +  geom_segment(aes(x = institution, y = preffalse.ind.l95ci, xend = institution, yend = preffalse.ind.u95ci), color="grey60", alpha=.5, data = data.china.w6.institution) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + geom_hline(aes(yintercept=0.00), colour="grey40", linetype="dotted") 
data.china.w6.institution$variable<-"Institution"
colnames(data.china.w6.institution)[1] <- "level"

data.china.w6.class<-aggregate(wvs.china.w6$conf.gov.ob, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)
data.china.w6.class[,3]<-aggregate(wvs.china.w6$conf.gov.dk, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,4]<-aggregate(wvs.china.w6$resp.hr.ob, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,5]<-aggregate(wvs.china.w6$resp.hr.dk, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,6]<-aggregate(wvs.china.w6$democraticness.ob, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,7]<-aggregate(wvs.china.w6$democraticness.dk, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,8]<-aggregate(wvs.china.w6$conf.tv.ob, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,9]<-aggregate(wvs.china.w6$conf.tv.dk, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,10]<-aggregate(wvs.china.w6$gentrust.ob, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,11]<-aggregate(wvs.china.w6$gentrust.dk, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,12]<-aggregate(wvs.china.w6$lifesat.ob, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,13]<-aggregate(wvs.china.w6$lifesat.dk, by=list(wvs.china.w6$class), FUN="sum", na.rm=TRUE)[2]
data.china.w6.class[,14]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$class), FUN="mean", na.rm=TRUE)[2]
data.china.w6.class[,15]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$class), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w6.class)<-c("class","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w6.class$conf.gov.dk.rate<-data.china.w6.class$conf.gov.dk/data.china.w6.class$conf.gov.ob
data.china.w6.class$resp.hr.dk.rate<-data.china.w6.class$resp.hr.dk/data.china.w6.class$resp.hr.ob
data.china.w6.class$democraticness.dk.rate<-data.china.w6.class$democraticness.dk/data.china.w6.class$democraticness.ob
data.china.w6.class$conf.tv.dk.rate<-data.china.w6.class$conf.tv.dk/data.china.w6.class$conf.tv.ob
data.china.w6.class$gentrust.dk.rate<-data.china.w6.class$gentrust.dk/data.china.w6.class$gentrust.ob
data.china.w6.class$lifesat.dk.rate<-data.china.w6.class$lifesat.dk/data.china.w6.class$lifesat.ob
data.china.w6.class$sens.ob<-data.china.w6.class$conf.gov.ob+data.china.w6.class$resp.hr.ob+data.china.w6.class$democraticness.ob
data.china.w6.class$sens.ob[data.china.w6.class$conf.gov.ob==0]<-NA
data.china.w6.class$sens.ob[data.china.w6.class$resp.hr.ob==0]<-NA
data.china.w6.class$sens.ob[data.china.w6.class$democraticness.ob==0]<-NA
data.china.w6.class$sens.dk<-data.china.w6.class$conf.gov.dk+data.china.w6.class$resp.hr.dk+data.china.w6.class$democraticness.dk
data.china.w6.class$sens.dk[data.china.w6.class$conf.gov.ob==0]<-NA
data.china.w6.class$sens.dk[data.china.w6.class$resp.hr.ob==0]<-NA
data.china.w6.class$sens.dk[data.china.w6.class$democraticness.ob==0]<-NA
data.china.w6.class$sens.dk.rate<-data.china.w6.class$sens.dk/data.china.w6.class$sens.ob
data.china.w6.class$sens.dk.se<-sqrt((data.china.w6.class$sens.dk.rate*(1-data.china.w6.class$sens.dk.rate))/data.china.w6.class$sens.ob)
data.china.w6.class$sens.dk.u95ci<-data.china.w6.class$sens.dk.rate+1.96*data.china.w6.class$sens.dk.se
data.china.w6.class$sens.dk.l95ci<-data.china.w6.class$sens.dk.rate-1.96*data.china.w6.class$sens.dk.se
data.china.w6.class$nonsens.ob<-data.china.w6.class$conf.tv.ob+data.china.w6.class$lifesat.ob+data.china.w6.class$gentrust.ob
data.china.w6.class$nonsens.ob[data.china.w6.class$conf.tv.ob==0]<-NA
data.china.w6.class$nonsens.ob[data.china.w6.class$lifesat.ob==0]<-NA
data.china.w6.class$nonsens.ob[data.china.w6.class$gentrust.ob==0]<-NA
data.china.w6.class$nonsens.dk<-data.china.w6.class$conf.tv.dk+data.china.w6.class$lifesat.dk+data.china.w6.class$gentrust.dk
data.china.w6.class$nonsens.dk[data.china.w6.class$conf.tv.ob==0]<-NA
data.china.w6.class$nonsens.dk[data.china.w6.class$lifesat.ob==0]<-NA
data.china.w6.class$nonsens.dk[data.china.w6.class$gentrust.ob==0]<-NA
data.china.w6.class$nonsens.dk.rate<-data.china.w6.class$nonsens.dk/data.china.w6.class$nonsens.ob
data.china.w6.class$nonsens.dk.se<-sqrt((data.china.w6.class$nonsens.dk.rate*(1-data.china.w6.class$nonsens.dk.rate))/data.china.w6.class$nonsens.ob)
data.china.w6.class$nonsens.dk.u95ci<-data.china.w6.class$nonsens.dk.rate+1.96*data.china.w6.class$nonsens.dk.se
data.china.w6.class$nonsens.dk.l95ci<-data.china.w6.class$nonsens.dk.rate-1.96*data.china.w6.class$nonsens.dk.se
data.china.w6.class$preffalse.ind<-data.china.w6.class$sens.dk.rate-data.china.w6.class$nonsens.dk.rate
data.china.w6.class$preffalse.ind.se<-sqrt((data.china.w6.class$nonsens.dk.se^2)+(data.china.w6.class$sens.dk.se^2))
data.china.w6.class$preffalse.ind.u95ci<-data.china.w6.class$preffalse.ind+1.96*data.china.w6.class$preffalse.ind.se
data.china.w6.class$preffalse.ind.l95ci<-data.china.w6.class$preffalse.ind-1.96*data.china.w6.class$preffalse.ind.se
data.china.w6.class$conf.gov.se<-data.china.w6.class$conf.gov.sd/sqrt(data.china.w6.class$conf.gov.ob)
data.china.w6.class$conf.gov.u95ci<-data.china.w6.class$conf.gov.mean+1.96*data.china.w6.class$conf.gov.se
data.china.w6.class$conf.gov.l95ci<-data.china.w6.class$conf.gov.mean-1.96*data.china.w6.class$conf.gov.se
ggplot(data.china.w6.class, aes(x=class, y=preffalse.ind)) + xlab("class") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = FALSE)+ theme_bw() + theme(legend.title=element_blank()) + geom_point(alpha=.5)  +  geom_segment(aes(x = class, y = preffalse.ind.l95ci, xend = class, yend = preffalse.ind.u95ci), color="grey60", alpha=.5, data = data.china.w6.class) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + geom_hline(aes(yintercept=0.00), colour="grey40", linetype="dotted") 
data.china.w6.class$variable<-"Class"
colnames(data.china.w6.class)[1] <- "level"

data.china.w6<-rbind(data.china.w6.education, data.china.w6.gender, data.china.w6.age, data.china.w6.class, data.china.w6.institution)

data.china.w5$wave<-"WVS - Wave 5 (2007)"
data.china.w6$wave<-"WVS - Wave 6 (2012)"
data.china.w5w6<-rbind(data.china.w5, data.china.w6)

#China Survey#

load("In Search of Preference Falsification - China Survey 2008 Raw.Rdata")

###CODE SENSITIVE QUESTIONS###

data.cs$democracy.obs<-1
data.cs$democracy.dk<-NA
data.cs$democracy.dk[is.na(data.cs$c1e)==TRUE]<-1
data.cs$democracy.dk[data.cs$c1e=="1"]<-0
data.cs$democracy.dk[data.cs$c1e=="2"]<-0
data.cs$democracy.dk[data.cs$c1e=="3"]<-0
data.cs$democracy.dk[data.cs$c1e=="4"]<-0
data.cs$democracy.dk[data.cs$c1e=="5"]<-0
data.cs$democracy.dk[data.cs$c1e=="6"]<-0
data.cs$democracy.dk[data.cs$c1e=="7"]<-0
data.cs$democracy.dk[data.cs$c1e=="8"]<-0
data.cs$democracy.dk[data.cs$c1e=="9"]<-0
data.cs$democracy.dk[data.cs$c1e=="10"]<-0
ftable(data.cs$democracy.dk)
ftable(data.cs$democracy.obs)

data.cs$humanrights.obs<-1
data.cs$humanrights.dk<-NA
data.cs$humanrights.dk[is.na(data.cs$d31g)==TRUE]<-1
data.cs$humanrights.dk[data.cs$d31g=="dont know"]<-1
data.cs$humanrights.dk[data.cs$d31g=="no answer"]<-1
data.cs$humanrights.dk[data.cs$d31g=="agree strongly"]<-0
data.cs$humanrights.dk[data.cs$d31g=="agree"]<-0
data.cs$humanrights.dk[data.cs$d31g=="neither agree nor disagree"]<-0
data.cs$humanrights.dk[data.cs$d31g=="disagree"]<-0
data.cs$humanrights.dk[data.cs$d31g=="disagree strongly"]<-0
ftable(data.cs$humanrights.dk)
ftable(data.cs$humanrights.obs)

data.cs$centgovsat.obs<-1
data.cs$centgovsat.dk<-NA
data.cs$centgovsat.dk[is.na(data.cs$h2b)==TRUE]<-1
data.cs$centgovsat.dk[data.cs$c9a=="0"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="1"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="2"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="3"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="4"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="5"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="6"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="7"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="8"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="9"]<-0
data.cs$centgovsat.dk[data.cs$c9a=="10"]<-0
ftable(data.cs$centgovsat.dk)
ftable(data.cs$centgovsat.obs)

data.cs$gentrust.obs<-1
data.cs$gentrust.dk<-NA
data.cs$gentrust.dk[data.cs$b6=="don't know"]<-1
data.cs$gentrust.dk[data.cs$b6=="no answer"]<-1
data.cs$gentrust.dk[is.na(data.cs$b6)==TRUE]<-1
data.cs$gentrust.dk[data.cs$b6=="people can be trusted"]<-0
data.cs$gentrust.dk[data.cs$b6=="you can't be too careful"]<-0
ftable(data.cs$gentrust.dk)
ftable(data.cs$gentrust.obs)

data.cs$lifesat.obs<-1
data.cs$lifesat.dk<-NA
data.cs$lifesat.dk[is.na(data.cs$h2b)==TRUE]<-1
data.cs$lifesat.dk[data.cs$h2b=="1"]<-0
data.cs$lifesat.dk[data.cs$h2b=="2"]<-0
data.cs$lifesat.dk[data.cs$h2b=="3"]<-0
data.cs$lifesat.dk[data.cs$h2b=="4"]<-0
data.cs$lifesat.dk[data.cs$h2b=="5"]<-0
data.cs$lifesat.dk[data.cs$h2b=="6"]<-0
data.cs$lifesat.dk[data.cs$h2b=="7"]<-0
data.cs$lifesat.dk[data.cs$h2b=="8"]<-0
data.cs$lifesat.dk[data.cs$h2b=="9"]<-0
data.cs$lifesat.dk[data.cs$h2b=="10"]<-0
ftable(data.cs$lifesat.dk)
ftable(data.cs$lifesat.obs)

data.cs$education.obs<-1
data.cs$education.dk<-NA
data.cs$education.dk[is.na(data.cs$c1b)==TRUE]<-1
data.cs$education.dk[data.cs$c1b=="0"]<-0
data.cs$education.dk[data.cs$c1b=="1"]<-0
data.cs$education.dk[data.cs$c1b=="2"]<-0
data.cs$education.dk[data.cs$c1b=="3"]<-0
data.cs$education.dk[data.cs$c1b=="4"]<-0
data.cs$education.dk[data.cs$c1b=="5"]<-0
data.cs$education.dk[data.cs$c1b=="6"]<-0
data.cs$education.dk[data.cs$c1b=="7"]<-0
data.cs$education.dk[data.cs$c1b=="8"]<-0
data.cs$education.dk[data.cs$c1b=="9"]<-0
data.cs$education.dk[data.cs$c1b=="10"]<-0
ftable(data.cs$education.dk)
ftable(data.cs$education.obs)

###CREATE AGGREGATE DATASETS###

data.cs$birthyear2<-2008-data.cs$age

data.cs$birthyear<-NA
data.cs$birthyear[data.cs$birthyear2<2000]<-"1980+"
data.cs$birthyear[data.cs$birthyear2<1980]<-"1970s"
data.cs$birthyear[data.cs$birthyear2<1970]<-"1960s"
data.cs$birthyear[data.cs$birthyear2<1960]<-"1950s"
data.cs$birthyear[data.cs$birthyear2<1950]<-"<1950"
ftable(data.cs$birthyear)

data.cs$class<-NA
data.cs$class[data.cs$z6=="lower class"]<-"1 - Lower class"
#data.cs$class[data.cs$z6=="Working class"]<-"2 - Lower middle class"
data.cs$class[data.cs$z6=="middle class"]<-"3 - Middle class"
data.cs$class[data.cs$z6=="upper-middle class"]<-"4 - Upper middle class"
data.cs$class[data.cs$z6=="upper class"]<-"5 - Upper class"
ftable(data.cs$class)

data.cs$education<-NA
data.cs$education<-"No education"
data.cs$education[data.cs$f16a>0]<-"Elementary"
data.cs$education[data.cs$f16b>0]<-"Secondary"
data.cs$education[data.cs$f16c>0]<-"Secondary"
data.cs$education[data.cs$f16d>0]<-"Secondary"
data.cs$education[data.cs$f16e>0]<-"University"
data.cs$education[data.cs$f16f>0]<-"University"
data.cs$education[data.cs$f16g>0]<-"University"
data.cs$education[data.cs$f16h>0]<-"University"
ftable(data.cs$education)

data.cs$institution<-"Private"
data.cs$institution[data.cs$f18=="employee of government agency, party agency, social organiza"]<-"Public"
data.cs$institution[data.cs$f18=="serviceman or police officer"]<-"Public"
data.cs$institution[data.cs$f18=="no answer"]<-NA
data.cs$institution[data.cs$f18=="do not konw"]<-NA #this typo is contained in the survey
data.cs$institution[data.cs$f20=="state-owned firm"]<-"Public"
data.cs$institution[data.cs$f20=="mostly state-owned enterprises"]<-"Public"
data.cs$institution[data.cs$f20=="party and government organizations"]<-"Public"
data.cs$institution[data.cs$f20=="state-owned institutional unit"]<-"Public"
data.cs$institution[data.cs$f20=="individually-owned business"]<-"Private"
data.cs$institution[data.cs$f20=="town or village owned enterprise"]<-"Public"
data.cs$institution[data.cs$f20=="private-owned firm"]<-"Private"
data.cs$institution[data.cs$f20=="institutional unit"]<-"Private"
data.cs$institution[data.cs$f20=="state-owned institutional unit"]<-"Public"
data.cs$institution[data.cs$f18=="no answer"]<-NA
data.cs$institution[data.cs$f18=="don't know"]<-NA #this typo is contained in the survey
ftable(data.cs$institution)

data.cs$female<-"Male"
data.cs$female[data.cs$gender=="female"]<-"Female"

data.cs$minority<-NA
data.cs$minority[data.cs$a8=="han"]<-"Han"
data.cs$minority[data.cs$a8=="other"]<-"Minority"

data.cs$rural<-NA
data.cs$rural[data.cs$a9=="urban"]<-"Urban"
data.cs$rural[data.cs$a9=="rural"]<-"Rural"

data.cs$ccp<-NA
data.cs$ccp[data.cs$e6_a=="yes"]<-"CCP Member"
data.cs$ccp[data.cs$e6_a=="no"]<-"Non Member"
data.cs$ccp[data.cs$e6_a=="don't know"]<-NA
data.cs$ccp[data.cs$e6_a=="no answer"]<-NA
data.cs$ccp[data.cs$e6_a=="no participation"]<-"Non Member"
data.cs$ccp[data.cs$e6_a=="once was"]<-"CCP Member"

vars<-c("birthyear", "birthyear2","class","education","institution","female", "minority", "rural", "ccp")

for (k in vars) {
  eval(parse(text=paste("data.cs.",k,"<-aggregate(data.cs$democracy.ob, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,3]<-aggregate(data.cs$democracy.dk, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,4]<-aggregate(data.cs$humanrights.ob, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,5]<-aggregate(data.cs$humanrights.dk, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,6]<-aggregate(data.cs$centgovsat.ob, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,7]<-aggregate(data.cs$centgovsat.dk, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,8]<-aggregate(data.cs$gentrust.ob, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,9]<-aggregate(data.cs$gentrust.dk, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,10]<-aggregate(data.cs$lifesat.ob, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,11]<-aggregate(data.cs$lifesat.dk, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,12]<-aggregate(data.cs$education.ob, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,13]<-aggregate(data.cs$education.dk, by=list(data.cs$",k,"), FUN='sum', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,14]<-aggregate(data.cs$c9a, by=list(data.cs$",k,"), FUN='mean', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("data.cs.",k,"[,15]<-aggregate(data.cs$c9a, by=list(data.cs$",k,"), FUN='sd', na.rm=TRUE)[2]",sep="")))
  eval(parse(text=paste("colnames(data.cs.",k,")<-c('label','democracy.ob','democracy.dk','humanrights.ob','humanrights.dk','centgovsat.ob','centgovsat.dk','gentrust.ob','gentrust.dk','lifesat.ob','lifesat.dk','education.ob','education.dk','centgovsat.mean','centgovsat.sd')",sep="")))
  eval(parse(text=paste("data.cs.",k,"$democracy.dk.rate<-data.cs.",k,"$democracy.dk/data.cs.",k,"$democracy.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$humanrights.dk.rate<-data.cs.",k,"$humanrights.dk/data.cs.",k,"$humanrights.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$centgovsat.dk.rate<-data.cs.",k,"$centgovsat.dk/data.cs.",k,"$centgovsat.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$education.dk.rate<-data.cs.",k,"$education.dk/data.cs.",k,"$education.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$gentrust.dk.rate<-data.cs.",k,"$gentrust.dk/data.cs.",k,"$gentrust.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$lifesat.dk.rate<-data.cs.",k,"$lifesat.dk/data.cs.",k,"$lifesat.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$sens.ob<-data.cs.",k,"$democracy.ob+data.cs.",k,"$humanrights.ob+data.cs.",k,"$centgovsat.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$sens.dk<-data.cs.",k,"$democracy.dk+data.cs.",k,"$humanrights.dk+data.cs.",k,"$centgovsat.dk",sep="")))
  eval(parse(text=paste("data.cs.",k,"$sens.dk.rate<-data.cs.",k,"$sens.dk/data.cs.",k,"$sens.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$sens.dk.se<-sqrt((data.cs.",k,"$sens.dk.rate*(1-data.cs.",k,"$sens.dk.rate))/data.cs.",k,"$sens.ob)",sep="")))
  eval(parse(text=paste("data.cs.",k,"$sens.dk.u95ci<-data.cs.",k,"$sens.dk.rate+1.96*data.cs.",k,"$sens.dk.se",sep="")))
  eval(parse(text=paste("data.cs.",k,"$sens.dk.l95ci<-data.cs.",k,"$sens.dk.rate-1.96*data.cs.",k,"$sens.dk.se",sep="")))
  eval(parse(text=paste("data.cs.",k,"$nonsens.ob<-data.cs.",k,"$education.ob+data.cs.",k,"$lifesat.ob+data.cs.",k,"$gentrust.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$nonsens.dk<-data.cs.",k,"$education.dk+data.cs.",k,"$lifesat.dk+data.cs.",k,"$gentrust.dk",sep="")))
  eval(parse(text=paste("data.cs.",k,"$nonsens.dk.rate<-data.cs.",k,"$nonsens.dk/data.cs.",k,"$nonsens.ob",sep="")))
  eval(parse(text=paste("data.cs.",k,"$nonsens.dk.se<-sqrt((data.cs.",k,"$nonsens.dk.rate*(1-data.cs.",k,"$nonsens.dk.rate))/data.cs.",k,"$nonsens.ob)",sep="")))
  eval(parse(text=paste("data.cs.",k,"$nonsens.dk.u95ci<-data.cs.",k,"$nonsens.dk.rate+1.96*data.cs.",k,"$nonsens.dk.se",sep="")))
  eval(parse(text=paste("data.cs.",k,"$nonsens.dk.l95ci<-data.cs.",k,"$nonsens.dk.rate-1.96*data.cs.",k,"$nonsens.dk.se",sep="")))
  eval(parse(text=paste("data.cs.",k,"$preffalse.ind<-data.cs.",k,"$sens.dk.rate-data.cs.",k,"$nonsens.dk.rate",sep="")))
  eval(parse(text=paste("data.cs.",k,"$preffalse.ind.se<-sqrt((data.cs.",k,"$nonsens.dk.se^2)+(data.cs.",k,"$sens.dk.se^2))",sep="")))
  eval(parse(text=paste("data.cs.",k,"$preffalse.ind.u95ci<-data.cs.",k,"$preffalse.ind+1.96*data.cs.",k,"$preffalse.ind.se",sep="")))
  eval(parse(text=paste("data.cs.",k,"$preffalse.ind.l95ci<-data.cs.",k,"$preffalse.ind-1.96*data.cs.",k,"$preffalse.ind.se",sep="")))
  eval(parse(text=paste("data.cs.",k,"$centgovsat.se<-data.cs.",k,"$centgovsat.sd/sqrt(data.cs.",k,"$centgovsat.ob)",sep="")))
  eval(parse(text=paste("data.cs.",k,"$centgovsat.u95ci<-data.cs.",k,"$centgovsat.mean+1.96*data.cs.",k,"$centgovsat.se",sep="")))
  eval(parse(text=paste("data.cs.",k,"$centgovsat.l95ci<-data.cs.",k,"$centgovsat.mean-1.96*data.cs.",k,"$centgovsat.se",sep="")))
  eval(parse(text=paste("data.cs.",k,"$variable<-'",k,"'",sep="")))
}

data.cs.var<-rbind(data.cs.birthyear,data.cs.class,data.cs.education,data.cs.institution,data.cs.ccp,data.cs.female,data.cs.minority,data.cs.rural)

data.cs.var$variable[data.cs.var$variable=="minority"]<-"Ethnicity"
data.cs.var$variable[data.cs.var$variable=="birthyear"]<-"Birthyear"
data.cs.var$variable[data.cs.var$variable=="ccp"]<-"Party"
data.cs.var$variable[data.cs.var$variable=="female"]<-"Gender"
data.cs.var$variable[data.cs.var$variable=="rural"]<-"Hukou"
data.cs.var$variable[data.cs.var$variable=="class"]<-"Class"
data.cs.var$variable[data.cs.var$variable=="education"]<-"Education"
data.cs.var$variable[data.cs.var$variable=="institution"]<-"Institution"

###MERGE IN WITH WVS AGGREGATES###

cs.figure<-data.frame(cbind(data.cs.var$label, data.cs.var$preffalse.ind, data.cs.var$preffalse.ind.l95ci, data.cs.var$preffalse.ind.u95ci, data.cs.var$variable))
colnames(cs.figure)<-c("level","preffalse.ind","preffalse.ind.l95ci","preffalse.ind.u95ci","variable")
cs.figure$wave<-"China Survey (2008)"

figure<-data.frame(cbind(data.china.w5w6$level, data.china.w5w6$preffalse.ind, data.china.w5w6$preffalse.ind.l95ci, data.china.w5w6$preffalse.ind.u95ci, data.china.w5w6$variable, data.china.w5w6$wave))
colnames(figure)<-c("level","preffalse.ind","preffalse.ind.l95ci","preffalse.ind.u95ci","variable","wave")

figure<-rbind(figure,cs.figure)
figure$preffalse.ind<-as.numeric(paste(figure$preffalse.ind))
figure$preffalse.ind.l95ci<-as.numeric(paste(figure$preffalse.ind.l95ci))
figure$preffalse.ind.u95ci<-as.numeric(paste(figure$preffalse.ind.u95ci))

###FIGURES###

pdf('fig5-china-agg.pdf', width=7.5, height=5.5)                                                                                                                                                                                                                                                                                                                                                                                                                                              
ggplot(figure, aes(y=level, x=preffalse.ind)) + ylab("Variable") + xlab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.8, face="bold")) + scale_shape(solid = TRUE)+ theme_classic() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue4", size=1.75)  +  geom_segment(aes(y = level, x = preffalse.ind.l95ci, yend = level, xend = preffalse.ind.u95ci), color="dodgerblue4", alpha=.4, lwd=.7, data = figure) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + xlim(c(-0.05,.30))+facet_grid(variable ~ wave, scales = "free", space = "free") + theme(strip.text.y = element_text(angle = 0))  +  theme_minimal() +theme(strip.text.y = element_text(angle = 360)) + geom_vline(xintercept = -0.05, size=.5,  color="grey20", alpha=.3) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))
dev.off()
                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         dev.off()

data.china.w5.age<-aggregate(wvs.china.w5$conf.gov.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)
data.china.w5.age[,3]<-aggregate(wvs.china.w5$conf.gov.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,4]<-aggregate(wvs.china.w5$resp.hr.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,5]<-aggregate(wvs.china.w5$resp.hr.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,6]<-aggregate(wvs.china.w5$democraticness.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,7]<-aggregate(wvs.china.w5$democraticness.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,8]<-aggregate(wvs.china.w5$conf.tv.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,9]<-aggregate(wvs.china.w5$conf.tv.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,10]<-aggregate(wvs.china.w5$gentrust.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,11]<-aggregate(wvs.china.w5$gentrust.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,12]<-aggregate(wvs.china.w5$lifesat.ob, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,13]<-aggregate(wvs.china.w5$lifesat.dk, by=list(wvs.china.w5$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w5.age[,14]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$birthyear), FUN="mean", na.rm=TRUE)[2]
data.china.w5.age[,15]<-aggregate(wvs.china.w5$conf.gov, by=list(wvs.china.w5$birthyear), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w5.age)<-c("birthyear","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w5.age$conf.gov.dk.rate<-data.china.w5.age$conf.gov.dk/data.china.w5.age$conf.gov.ob
data.china.w5.age$resp.hr.dk.rate<-data.china.w5.age$resp.hr.dk/data.china.w5.age$resp.hr.ob
data.china.w5.age$democraticness.dk.rate<-data.china.w5.age$democraticness.dk/data.china.w5.age$democraticness.ob
data.china.w5.age$conf.tv.dk.rate<-data.china.w5.age$conf.tv.dk/data.china.w5.age$conf.tv.ob
data.china.w5.age$gentrust.dk.rate<-data.china.w5.age$gentrust.dk/data.china.w5.age$gentrust.ob
data.china.w5.age$lifesat.dk.rate<-data.china.w5.age$lifesat.dk/data.china.w5.age$lifesat.ob
data.china.w5.age$sens.ob<-data.china.w5.age$conf.gov.ob+data.china.w5.age$resp.hr.ob+data.china.w5.age$democraticness.ob
data.china.w5.age$sens.ob[data.china.w5.age$conf.gov.ob==0]<-NA
data.china.w5.age$sens.ob[data.china.w5.age$resp.hr.ob==0]<-NA
data.china.w5.age$sens.ob[data.china.w5.age$democraticness.ob==0]<-NA
data.china.w5.age$sens.dk<-data.china.w5.age$conf.gov.dk+data.china.w5.age$resp.hr.dk+data.china.w5.age$democraticness.dk
data.china.w5.age$sens.dk[data.china.w5.age$conf.gov.ob==0]<-NA
data.china.w5.age$sens.dk[data.china.w5.age$resp.hr.ob==0]<-NA
data.china.w5.age$sens.dk[data.china.w5.age$democraticness.ob==0]<-NA
data.china.w5.age$sens.dk.rate<-data.china.w5.age$sens.dk/data.china.w5.age$sens.ob
data.china.w5.age$sens.dk.se<-sqrt((data.china.w5.age$sens.dk.rate*(1-data.china.w5.age$sens.dk.rate))/data.china.w5.age$sens.ob)
data.china.w5.age$sens.dk.u95ci<-data.china.w5.age$sens.dk.rate+1.96*data.china.w5.age$sens.dk.se
data.china.w5.age$sens.dk.l95ci<-data.china.w5.age$sens.dk.rate-1.96*data.china.w5.age$sens.dk.se
data.china.w5.age$nonsens.ob<-data.china.w5.age$conf.tv.ob+data.china.w5.age$lifesat.ob+data.china.w5.age$gentrust.ob
data.china.w5.age$nonsens.ob[data.china.w5.age$conf.tv.ob==0]<-NA
data.china.w5.age$nonsens.ob[data.china.w5.age$lifesat.ob==0]<-NA
data.china.w5.age$nonsens.ob[data.china.w5.age$gentrust.ob==0]<-NA
data.china.w5.age$nonsens.dk<-data.china.w5.age$conf.tv.dk+data.china.w5.age$lifesat.dk+data.china.w5.age$gentrust.dk
data.china.w5.age$nonsens.dk[data.china.w5.age$conf.tv.ob==0]<-NA
data.china.w5.age$nonsens.dk[data.china.w5.age$lifesat.ob==0]<-NA
data.china.w5.age$nonsens.dk[data.china.w5.age$gentrust.ob==0]<-NA
data.china.w5.age$nonsens.dk.rate<-data.china.w5.age$nonsens.dk/data.china.w5.age$nonsens.ob
data.china.w5.age$nonsens.dk.se<-sqrt((data.china.w5.age$nonsens.dk.rate*(1-data.china.w5.age$nonsens.dk.rate))/data.china.w5.age$nonsens.ob)
data.china.w5.age$nonsens.dk.u95ci<-data.china.w5.age$nonsens.dk.rate+1.96*data.china.w5.age$nonsens.dk.se
data.china.w5.age$nonsens.dk.l95ci<-data.china.w5.age$nonsens.dk.rate-1.96*data.china.w5.age$nonsens.dk.se
data.china.w5.age$preffalse.ind<-data.china.w5.age$sens.dk.rate-data.china.w5.age$nonsens.dk.rate
data.china.w5.age$preffalse.ind.se<-sqrt((data.china.w5.age$nonsens.dk.se^2)+(data.china.w5.age$sens.dk.se^2))
data.china.w5.age$preffalse.ind.u95ci<-data.china.w5.age$preffalse.ind+1.96*data.china.w5.age$preffalse.ind.se
data.china.w5.age$preffalse.ind.l95ci<-data.china.w5.age$preffalse.ind-1.96*data.china.w5.age$preffalse.ind.se
data.china.w5.age$conf.gov.se<-data.china.w5.age$conf.gov.sd/sqrt(data.china.w5.age$conf.gov.ob)
data.china.w5.age$conf.gov.u95ci<-data.china.w5.age$conf.gov.mean+1.96*data.china.w5.age$conf.gov.se
data.china.w5.age$conf.gov.l95ci<-data.china.w5.age$conf.gov.mean-1.96*data.china.w5.age$conf.gov.se
data.china.w5.age$variable<-"Birthyear"
colnames(data.china.w5.age)[1] <- "level"

data.china.w6.age<-aggregate(wvs.china.w6$conf.gov.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)
data.china.w6.age[,3]<-aggregate(wvs.china.w6$conf.gov.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,4]<-aggregate(wvs.china.w6$resp.hr.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,5]<-aggregate(wvs.china.w6$resp.hr.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,6]<-aggregate(wvs.china.w6$democraticness.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,7]<-aggregate(wvs.china.w6$democraticness.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,8]<-aggregate(wvs.china.w6$conf.tv.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,9]<-aggregate(wvs.china.w6$conf.tv.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,10]<-aggregate(wvs.china.w6$gentrust.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,11]<-aggregate(wvs.china.w6$gentrust.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,12]<-aggregate(wvs.china.w6$lifesat.ob, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,13]<-aggregate(wvs.china.w6$lifesat.dk, by=list(wvs.china.w6$birthyear), FUN="sum", na.rm=TRUE)[2]
data.china.w6.age[,14]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$birthyear), FUN="mean", na.rm=TRUE)[2]
data.china.w6.age[,15]<-aggregate(wvs.china.w6$conf.gov, by=list(wvs.china.w6$birthyear), FUN="sd", na.rm=TRUE)[2]
colnames(data.china.w6.age)<-c("birthyear","conf.gov.ob","conf.gov.dk","resp.hr.ob","resp.hr.dk","democraticness.ob","democraticness.dk","conf.tv.ob","conf.tv.dk","gentrust.ob","gentrust.dk","lifesat.ob","lifesat.dk","conf.gov.mean","conf.gov.sd")
data.china.w6.age$conf.gov.dk.rate<-data.china.w6.age$conf.gov.dk/data.china.w6.age$conf.gov.ob
data.china.w6.age$resp.hr.dk.rate<-data.china.w6.age$resp.hr.dk/data.china.w6.age$resp.hr.ob
data.china.w6.age$democraticness.dk.rate<-data.china.w6.age$democraticness.dk/data.china.w6.age$democraticness.ob
data.china.w6.age$conf.tv.dk.rate<-data.china.w6.age$conf.tv.dk/data.china.w6.age$conf.tv.ob
data.china.w6.age$gentrust.dk.rate<-data.china.w6.age$gentrust.dk/data.china.w6.age$gentrust.ob
data.china.w6.age$lifesat.dk.rate<-data.china.w6.age$lifesat.dk/data.china.w6.age$lifesat.ob
data.china.w6.age$sens.ob<-data.china.w6.age$conf.gov.ob+data.china.w6.age$resp.hr.ob+data.china.w6.age$democraticness.ob
data.china.w6.age$sens.ob[data.china.w6.age$conf.gov.ob==0]<-NA
data.china.w6.age$sens.ob[data.china.w6.age$resp.hr.ob==0]<-NA
data.china.w6.age$sens.ob[data.china.w6.age$democraticness.ob==0]<-NA
data.china.w6.age$sens.dk<-data.china.w6.age$conf.gov.dk+data.china.w6.age$resp.hr.dk+data.china.w6.age$democraticness.dk
data.china.w6.age$sens.dk[data.china.w6.age$conf.gov.ob==0]<-NA
data.china.w6.age$sens.dk[data.china.w6.age$resp.hr.ob==0]<-NA
data.china.w6.age$sens.dk[data.china.w6.age$democraticness.ob==0]<-NA
data.china.w6.age$sens.dk.rate<-data.china.w6.age$sens.dk/data.china.w6.age$sens.ob
data.china.w6.age$sens.dk.se<-sqrt((data.china.w6.age$sens.dk.rate*(1-data.china.w6.age$sens.dk.rate))/data.china.w6.age$sens.ob)
data.china.w6.age$sens.dk.u95ci<-data.china.w6.age$sens.dk.rate+1.96*data.china.w6.age$sens.dk.se
data.china.w6.age$sens.dk.l95ci<-data.china.w6.age$sens.dk.rate-1.96*data.china.w6.age$sens.dk.se
data.china.w6.age$nonsens.ob<-data.china.w6.age$conf.tv.ob+data.china.w6.age$lifesat.ob+data.china.w6.age$gentrust.ob
data.china.w6.age$nonsens.ob[data.china.w6.age$conf.tv.ob==0]<-NA
data.china.w6.age$nonsens.ob[data.china.w6.age$lifesat.ob==0]<-NA
data.china.w6.age$nonsens.ob[data.china.w6.age$gentrust.ob==0]<-NA
data.china.w6.age$nonsens.dk<-data.china.w6.age$conf.tv.dk+data.china.w6.age$lifesat.dk+data.china.w6.age$gentrust.dk
data.china.w6.age$nonsens.dk[data.china.w6.age$conf.tv.ob==0]<-NA
data.china.w6.age$nonsens.dk[data.china.w6.age$lifesat.ob==0]<-NA
data.china.w6.age$nonsens.dk[data.china.w6.age$gentrust.ob==0]<-NA
data.china.w6.age$nonsens.dk.rate<-data.china.w6.age$nonsens.dk/data.china.w6.age$nonsens.ob
data.china.w6.age$nonsens.dk.se<-sqrt((data.china.w6.age$nonsens.dk.rate*(1-data.china.w6.age$nonsens.dk.rate))/data.china.w6.age$nonsens.ob)
data.china.w6.age$nonsens.dk.u95ci<-data.china.w6.age$nonsens.dk.rate+1.96*data.china.w6.age$nonsens.dk.se
data.china.w6.age$nonsens.dk.l95ci<-data.china.w6.age$nonsens.dk.rate-1.96*data.china.w6.age$nonsens.dk.se
data.china.w6.age$preffalse.ind<-data.china.w6.age$sens.dk.rate-data.china.w6.age$nonsens.dk.rate
data.china.w6.age$preffalse.ind.se<-sqrt((data.china.w6.age$nonsens.dk.se^2)+(data.china.w6.age$sens.dk.se^2))
data.china.w6.age$preffalse.ind.u95ci<-data.china.w6.age$preffalse.ind+1.96*data.china.w6.age$preffalse.ind.se
data.china.w6.age$preffalse.ind.l95ci<-data.china.w6.age$preffalse.ind-1.96*data.china.w6.age$preffalse.ind.se
data.china.w6.age$conf.gov.se<-data.china.w6.age$conf.gov.sd/sqrt(data.china.w6.age$conf.gov.ob)
data.china.w6.age$conf.gov.u95ci<-data.china.w6.age$conf.gov.mean+1.96*data.china.w6.age$conf.gov.se
data.china.w6.age$conf.gov.l95ci<-data.china.w6.age$conf.gov.mean-1.96*data.china.w6.age$conf.gov.se
data.china.w6.age$variable<-"Birthyear"
colnames(data.china.w6.age)[1] <- "level"

pdf('fig5-china-age.pdf', width=6.5, height=6.5)  
data.china.w5.age$birthyear2<-as.numeric(paste(data.china.w5.age$level))
data.china.w6.age$birthyear2<-as.numeric(paste(data.china.w6.age$level))
p1<-ggplot(data.china.w5.age, aes(x=level, y=conf.gov.mean)) + xlab("Birthyear") + ylab("Confidence in Government (mean)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + scale_shape(solid = FALSE)+ theme_classic() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue4", size=1.75)  +  geom_segment(aes(x = level, y = conf.gov.l95ci, xend = level, yend = conf.gov.u95ci), color="dodgerblue4", alpha=.4, lwd=.7, data = data.china.w5.age) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + ggtitle("WVS - Wave 5 (2007)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) + theme(plot.title = element_text(size = 11))
p2<-ggplot(data.china.w5.age, aes(x=level, y=preffalse.ind)) + xlab("Birthyear") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + scale_shape(solid = FALSE)+ theme_classic() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue4", size=1.75)  +  geom_segment(aes(x = level, y = preffalse.ind.l95ci, xend = level, yend = preffalse.ind.u95ci), color="dodgerblue4", alpha=.4, lwd=.7, data = data.china.w5.age) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + ggtitle("WVS - Wave 5 (2007)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) + theme(plot.title = element_text(size = 11))  + ylim(c(-.0,.25))
p3<-ggplot(data.china.w6.age, aes(x=level, y=conf.gov.mean)) + xlab("Birthyear") + ylab("Confidence in Government (mean)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + scale_shape(solid = FALSE)+ theme_classic() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue4", size=1.75)  +  geom_segment(aes(x = level, y = conf.gov.l95ci, xend = level, yend = conf.gov.u95ci), color="dodgerblue4", alpha=.4, lwd=.7, data = data.china.w6.age) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + ggtitle("WVS - Wave 6 (2012)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) + theme(plot.title = element_text(size = 11))
p4<-ggplot(data.china.w6.age, aes(x=level, y=preffalse.ind)) + xlab("Birthyear") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + scale_shape(solid = FALSE)+ theme_classic() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue4", size=1.75)  +  geom_segment(aes(x = level, y = preffalse.ind.l95ci, xend = level, yend = preffalse.ind.u95ci), color="dodgerblue4", alpha=.4, lwd=.7, data = data.china.w6.age) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + ggtitle("WVS - Wave 6 (2012)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) + theme(plot.title = element_text(size = 11)) + ylim(c(-.0,.25))
p5<-ggplot(data.cs.birthyear, aes(x=label, y=centgovsat.mean)) + xlab("Birthyear") + ylab("Satisfaction with Government (mean)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + scale_shape(solid = FALSE)+ theme_classic() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue4", size=1.75)  +  geom_segment(aes(x = label, y = centgovsat.l95ci, xend = label, yend = centgovsat.u95ci), color="dodgerblue4", alpha=.4, lwd=.7, data = data.cs.birthyear) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + ggtitle("China Survey (2008)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) + theme(plot.title = element_text(size = 11)) + ylim(c(7,9))
p6<-ggplot(data.cs.birthyear, aes(x=label, y=preffalse.ind)) + xlab("Birthyear") + ylab("Self-censorship Index (mean)") + theme(plot.title = element_text(lineheight=.6, face="bold")) + scale_shape(solid = FALSE)+ theme_classic() + theme(legend.title=element_blank()) + geom_point(alpha=.5, color="dodgerblue4", size=1.75)  +  geom_segment(aes(x = label, y = preffalse.ind.l95ci, xend = label, yend = preffalse.ind.u95ci), color="dodgerblue4", alpha=.4, lwd=.7, data = data.cs.birthyear) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + geom_smooth(method=loess, se=FALSE, colour="black") + ggtitle("China Survey (2008)", subtitle = NULL) + theme(plot.title = element_text(hjust = 0.5)) + theme(plot.title = element_text(size = 11)) + ylim(c(-.0,.25))
arrange_ggplot2(p1,p3,p5,p2,p4,p6, nrow=2, ncol=3)
dev.off()










